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ABSTRACT 


The  goal  of  the  work  under  the  Conventional  Weapons  Effects  (cwe)  program  reported  here  is 
to  predict  the  dynamic  response  and  damage  of  buried  reinforced  concrete  structures  subjected 
to  airblast  and  bomb  case  fragment  impact  from  internal  detonations  of  conventional  weapons. 
Performing  such  calculations  with  a  first  principles  code  requires  implementation  of  a  comprehensive 
material  model  for  concrete,  particularly  for  modeling  strain-softening  in  the  tensile  regime  and  in 
the  low  confining  pressure  compressive  regime.  We  implemented  a  smooth-cap  model  with  damage 
into  the  nonlinear  finite  element  code  DYNA3D  (Reference  1)  to  analyze  the  dynamic  response  of 
plain  concrete,  and  other  geologic  materials.  The  capabilities  of  the  concrete  damage  model,  with 
and  without  modeling  steel  reinforcement,  are  demonstrated  in  the  present  report  through  a  series 
of  structural  benchmark  calculations.  The  validity  of  the  concrete  damage  and  reinforcement  models 
will  be  assessed  in  a  future  report  by  comparing  the  model  predictions  with  full  scale  Dipole  East 
tests  and  sub-scale  Precision  Model  tests. 

The  main  body  of  this  report  reviews  two  topics.  The  first  topic  is  concrete  material  model  develop¬ 
ment.  The  original  basis  of  our  concrete  damage  model  is  a  two-invariant,  smooth-cap  elasto-plastic 
model  (Reference  2).  We  added  numerous  featmes  to  this  model  to  improve  comparisons  with  stan¬ 
dard  laboratory  test  data.  These  improvements  include  a  damage  formulation  for  modeling  strain- 
softening  and  modulus  reduction  in  both  the  tensile  and  compression  regimes,  a  three-invariant 
plasticity  formulation  to  model  a  realistic  failure  envelope,  and  a  viscoplastic  formulation  for  mod¬ 
eling  strength  enhancement  at  high  strain  rates. 

The  second  topic  is  selected  benchmark  applications  using  the  concrete  damage  model.  This  re¬ 
port  includes  fits  of  the  model  to  standard  laboratory  test  data,  bending  analysis  of  reinforced 
concrete  slabs,  and  comparisons  of  single  and  multiple  element  simulations  for  unconfined  compres¬ 
sion  tests.  The  multi-element  test  simulations  predict  the  softening  response  and  localized  damage 
accumulation  of  concrete  test  specimens  better  them  the  single  element  calculations.  In  particular, 
multi-element  simulations  predict  more  severe  softening  than  single  element  models.  Multi-element 
simulations  also  predict  diagonal  damage  patterns  eind  splitting,  which  are  typical  failure,  models  of 
concrete  test  specimens. 

This  work  was  preformed  as  part  of  the  CWE  Program,  Loads  and  Dynamic  Response  of  Structures, 
under  contract  DNA001-91-C-0075.  The  DNA  technical  monitor  was  Mr.  Mike  Giltrud.  The  APTEK 
program  manager  was  Dr.  Brett  Lewis.  - 
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SECTION  1 

INTRODUCTION  TO  THE  CONCRETE  DAMAGE 

MODEL 


Numerous  plasticity  models  have  been  developed  over  the  years  to  describe  the  constitutive  behavior 
of  concrete  under  various  loading  conditions.  One  common  type  of  plasticity  model  is  the  elasto- 
plastic  cap  model  originally  attributed  to  Sandler  and  Rubin  in  Reference  3.  Cap  models  predict 
perfectly-plastic  response  for  laboratory  test  simulations,  such  as  triaxial  compression  (txc)  or 
extension  (txe)  simulations,  as  schematically  shown  by  the  line  with  circles  in  Figure  1.1.  Perfectly- 
plastic  response  is  typical  of  concrete  and  other  geologic  materials  at  high  confining  pressures,  but 
is  not  representative  of  concrete  at  low  or  no  confining  pressure.  At  low  confining  pressure,  concrete 
exhibits  strain-softening,  which  is  a  decrease  in  strength  during  progressive  straining  after  a  peak 
strength  value  is  reached. 


Axial  Strain 


Figure  1.1.  Comparison  of  strain-softening  (calculated  with  damage)  and  perfectly-plastic  (calcu¬ 
lated  without  damage)  responses  for  TXC  simulations. 

More  recently,  Simo  and  Ju  in  References  4-5  developed  a  continuum  damage  formulation  with 
application  to  cap  models,  which  captures  the  strain-softening  behavior  of  concrete.  APTEK  imple- 


1 


mented  their  damage  formulation  into  the  elasto-plastic  cap  model  in  the  DYNA3D  finite  element 
code.  A  simulation  with  strain-softening  is  indicated  by  the  solid  line  in  Figure  1.1.  The  cap  model 
with  damage  captures  the  essential  features  of  concrete  behavior,  such  cis  pre-peak  hardening,  post¬ 
peak  softening,  sheax-enhanced  compaction,  volume  expansion  under  compressive  loading  (dilation), 
modulus  reduction  under  cyclic  loading/unloading,  irreversible  deformation,  and  localized  damage 
accumulation,  as  discussed  throughout  the  remainder  of  this  report. 

The  cap  model  with  damage  treats  plastic  flow  and  damage  accumulation  as  separate  processes. 
Plastic  flow  is  generally  assumed  to  be  due  to  the  frictional  movement  (slip)  of  microcrack  surfaces 
and  is  controlled  by  the  presence  of  shear  stresses.  Plastic  flow  results  in  permanent  deforma¬ 
tion  although  the  elastic  moduli  do  not  degrade  with  this  modelling  mechanism.  Damage,  or  the 
progressive  degradation  of  moduli  and  strength  commonly  observed  in  concrete  test  specimens,  is 
thought  to  be  due  to  the  growth  and  coalescence  of  microcracks.  For  the  remainder  of  this  report, 
any  reference  to  damage  will  indicate  a  reduction  in  moduli  caused  by  the  growth  and  coalescence 
of  microcracks. 

The  significance  of  modeling  damage  accumulation  as  well  as  plastic  flow  is  best  explained  by 
examining  the  schematic  stress-strain  curves  shown  in  Figure  1.2  for  three  different  material  model 
behaviors.  The  elasto-plastic  model  in  Figure  1.2a  is  commonly  referred  to  as  a  plasticity  softening 
model.  The  stress-strain  simulation  exhibits  permanent  deformation  following  elastic  unloading, 
but  no  reduction  in  moduli.  Permanent  deformation  is  evident  because  the  simulation  unloads  to 
zero  stress  at  non-zero  values  of  strain.  Lack  of  damage  (modulus  reduction)  is  evident  because  the 
slopes  of  the  initial  loading  and  subsequent  loading/unloading  curves  axe  equal.  Such  behavior  is 
often  modeled  with  cap  models  by  shifting  or  contracting  the  failure  (plasticity)  surface  to  produce 
the  strain-softening  response.  APTEK,  however,  does  not  model  strain-softening  in  this  manner.  The 
second  model  in  Figure  1.2b  is  an  example  of  an  elasto-damage  model.  The  stress-strain  simulation 
exhibits  a  reduction  in  moduli  with  increasing  strain-softening,  but  no  permanent  deformation.  The 
reduction  in  moduli  is  indicated  by  the  decrease  in  the  elastic  loading/unloading  slopes  as  strain¬ 
softening  progresses.  Such  a  model  is  readily  implemented  by  applying  Simo  and  Ju’s  damage 
degradation  formulation  to  an  elastic  model,  rather  than  cap  model.  In  fact,  we  model  concrete 
response  in  the  brittle  (tensile  pressure)  regime  with  an  elasto-damage  simplification  of  the  cap 
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Figure  1.3.  Measured  cyclic  stress-strain  response  for  concrete  showing  degradation  of  the  elastic 
modulus  with  increasing  strain-softening. 

model.  The  third  model  in  Figure  1.2c  is  an  example  of  an  elasto-plastic  damage  model.  Both 
permanent  deformation  and  modulus  reduction  are  modeled.  APTEK’s  cap  model  with  damage 
described  in  this  report  is  an  elasto-plastic  damage  model,  and  is  used  for  modeling  concrete  response 
in  the  ductile  (compressive  pressure)  regime. 

Experimental  evidence  for  modulus  reduction  is  available  in  the  literature,  t.g.  the  cyclic  stress- 
strain  response  from  uniaxial  compressive  stress  tests  performed  on  concrete  by  Hurlbut  is  repro¬ 
duced  in  Figure  1.3  from  Reference  6.  Note  that  the  average  slope  from  the  loading/unloading 
loops  decreases  with  increasing  strain.  Experimental  evidence  for  modulus  reduction  from  lab¬ 
oratory  tests  conducted  by  WES  on  White  Sands  Missile  Range  (wsMR-5)  concrete  is  shown  in 
Figure  1.4.  In  these  tests,  concrete  cylinders  were  loaded  in  uniaxial  strain  to  a  pressure  of  138 
MPa  and  then  unloaded.  The  damaged  cylinder  was  then  loaded  in  either  direct  pull  or  unconfined 
compression.  Both  sets  of  data  indicate  a  substantial  reduction  in  loading  modulus  and  strength 
following  preloading  in  uniaxial  strain.  Also  shown  in  Figure  1.4  is  the  less  damaging  effect  of 
preloading  cylinders  in  isotropic  compression  (to  a  pressure  of  138  MPA)  prior  to  direct  pull  or 
unconfined  compression  tests. 

We  begin  our  discussion  of  the  cap  model  with  damage  with  a  brief  review  of  Simo  and  Ju’s  damage 
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Strain 


(a)  Elasto-plastic  or  plasticity  softening  model. 


(b)  Elasto-damage  model. 


(c)  Elasto-plastic  damage  model. 

Figure  1.2.  Comparison  of  material  model  behaviors. 
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(a)  Direct  pull  test  data. 


(b)  Unconfined  compression  test  data. 

Figure  1.4.  Reduction  in  modulus  and  strength  following  preloading  in  uniaxial  strain  for  WSMR-5 
concrete. 
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degradation  formulation  and  details  of  our  implementation  into  DYNA3D.  Implementations  for  mod¬ 
eling  ductile  and  brittle  response  are  discussed  separately  in  Sections  2  and  3.  Our  implementation 
for  ductile  response  combines  Simo  and  Ju’s  damage  degradation  formulation  with  the  smooth-cap 
model  implemented  into  DYNA3D  by  Pelessone  in  Reference  2,  a  three  stress-invariant  formulation 
developed  by  Rubin  in  Reference  7,  and  a  Duvaut-Lions  viscoplastic  update  algorithm  developed 
by  Simo  et.al.  in  Reference  12.  Our  implementation  for  brittle  response  is  a  three-invariant  elasto- 
damage  simplification  of  the  cap  model  in  which  plasticity  is  neglected.  Other  features  of  the 
implementation  include  modeling  bulk  modulus  ‘healing’  to  simulate  opening  and  closing  of  cracks, 
and  modeling  dilation  damage.  These  features  are  discussed  in  Sections  4  and  5,  respectively.  The 
implementation  of  these  formulations  into  a  cohesive  model  will  henceforth  be  referred  to  as  the 
APTEK  concrete  model. 

Verification  and  validation  of  the  concrete  model  is  in  progress  for  OWE  applications.  These  applica¬ 
tions  include  structural  benchmark  calculations  and  pre-test  predictions  for  Dipole  East  structures^. 
The  structural  benchmark  calculations  are  a  series  of  four  rounds  of  calculations,  performed  by  par¬ 
ticipating  analysts  on  the  CWE  program.  The  objectives  of  the  Round  A  and  B  calculations  are 
to  compare  and  assess  the  participants  models  for  plain  concrete.  These  calculations  examine  the 
effect  of  mesh  refinement  on  the  predicted  softening  response  and  accumulation  of  damage  for  TXC 
simulations.  The  objectives  of  the  Round  C  and  D  calculations  are  to  compare  the  participants 
models  for  reinforced  concrete.  These  calculations  study  the  dynamics  of  reinforced  slabs  sections. 
Comparisons  of  all  the  participants  computational  results  have  been  compiled  by  RDA  and  reported 
in  a  four  volumes  of  viewfoils  in  References  8  through  11.  Selected  benchmark  results  are  reported 
in  Section  6-8. 


‘Predictions  for  Dipole  East  structures  and  Round  E  benchmarks  will  be  reported  in  a  future  technicail  report. 
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SECTION  2 

MODELING  DAMAGE  IN  THE  DUCTILE  REGIME 


Simo  and  Ju  developed  alternative  strain-  and  stress-based  damage  models.  Damage  is  directly 
linked  to  the  history  of  strains  (stresses)  in  the  strain-based  (stress-based)  model.  The  strain-based 
and  stress-based  formulations  are  based  on  thermodynamic  considerations,  although  they  are  not 
equivalent  formulations. 

We  extended  the  elasto-plastic  cap  model  to  include  Simo  and  Ju’s  strain-based  damage  formulation 
for  modeling  ductile  response.  Our  terminology  for  ductile  response  is  any  state  of  stress  in  which 
the  pressure  is  compressive.  We  choose  to  implement  the  strain-based  formulation  rather  than  the 
stress-based  formulation  for  two  reasons.  First,  the  strain-based  formulation  is  relatively  easy  to 
implement,  as  discussed  in  subsequent  paragraphs,  because  it  does  not  involve  modification  of  the 
plasticity  return  mapping  algorithm.  Second,  we  felt  a  strain-based  damage  model  is  more  likely 
to  predict  a  “splitting  mode”,  i.e.  a  cylinder  subject  to  unconfined  uniaxial  compression  develops 
cracks  parallel  to  the  axis  of  loading,  due  to  the  presence  of  tensile  radial  and  hoop  strains. 

Although  only  the  strain-based  formulation  has  been  implemented  to  date,  future  efforts  should 
include  review,  and  possibly  implementation,  of  a  stress-based  damage  model  for  comparison  with 
the  strain-bcised  model.  The  relative  merits  and  drawbacks  of  each  type  of  formulation  are  of 
particular  interest  during  preloading/unloading  into  compression  (tension)  and  subsequent  loading 
into  tension  (compression). 

2.1  ISOTROPIC  FORMULATION. 

Simo  and  Ju  base  their  strain-based  damage  model  on  the  effective  stress  concept  and  the  hypothesis 
of  strain  equivalence.  This  hypothesis  states  that  “the  strain  aissociated  with  a  damaged  state  under 
the  applied  stress  is  equivalent  to  the  strain  associated  with  its  undamaged  state  under  the  effective 
stress.”  Let  a,-;  be  the  stress  tensor  associated  with  the  damaged  state,  and  a,-,-  be  the  effective 
stress  tensor,  or  stress  tensor  associated  with  the  undamaged  state.  A  scalar  damage  parameter  d 
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transforms  the  effective  stress  tensor  into  the  stress  tensor,  follows: 


cTij  =  (1  -  d)aij  (2.1) 

The  scalar  damage  parameter  d  ranges  from  zero  for  no  damage  and  approaches  unity  for  maximum 
damage.  Thus  1  —  d  is  a  reduction  factor  associated  with  the  amount  of  damage.  One  effect  of  the 
reduction  factor  1 —d  is  to  reduce  the  loading/unloading  slope,  or  moduli,  during  cyclic  loading.  For 
linear  elastic  response,  the  bulk  and  shecir  moduli,  K  and  G,  are  degraded  isotropically.  Another 
effect  is  to  model  the  softening  response  observed  in  lab  test  data.  A  physical  interpretation  of  the 
reduction  factor  1  —  d  is  the  ratio  of  undamaged  to  nominal  surface  area  at  a  material  point.  The 
effective  stress  d-,j  acts  on  the  undamaged  surface  area  while  the  total  stress  acts  on  the  nominal 
area. 

It  is  helpful  to  refer  back  to  the  unconfined  compression  simulations  previously  shown  in  Figure  1.1 
to  visualize  the  relationship  between  the  stress  associated  with  the  damaged  (strain-softening)  and 
undamaged  (perfectly-plastic)  states.  The  effective  stress  from  Equation  2.1  is  stress  calculated  for 
perfectly-plastic  response  in  Figure  1.1.  Thus  the  reduction  factor  1  —  d  decreases  with  axial  strain, 
and  is  equal  the  ratio  of  the  stress  calculated  with  damage  (strain-softening)  to  the  stress  calculated 
without  damage  (perfectly-plastic). 

2.2  EVOLUTION  OF  THE  DAMAGE  PARAMETER. 

It  is  assumed  that  the  damage  parameter  d  depends  on  the  history  of  total  strains.  In  particular,  d 
is  a  function  of  the  undamaged  energy  norm,  from  which  Simo  and  Ju  define  the  equivalent  strain, 
f,  as  follows: 


d  =  G{f) 
f  = 


(2.2) 

(2.3) 
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Here  is  the  undamaged  elastic  strain  energy.  For  the  linear  elastic  case,  =  \Cijki^ki^ij,  where 
Cijki  is  the  undamaged  linear  elasticity  tensor  and  are  the  total  strain  components.  One  way  of 
expanding  is: 

-  (o'lieii  +  0-22^22  +  <2‘|3e33)  +  <Ti2Cl2  +  +  <^23^23  (2.4) 

where  the  stress  tensor  is  calculated  from  the  total  strains  and  undamaged  moduli,  K  and  G. 

Damage  accumulates  when  f  exceeds  a  damage  threshold,  r,  i.e.  when  fn+i  >  r^.  Here  the  subscripts 
n  denote  the  time  step  in  the  finite  element  analysis.  The  initial  damage  threshold  is  denoted  tq. 
The  damage  threshold  at  each  time  step  is  held  equal  to  the  initial  damage  threshold  (rn  =  tq)  until 
fn+i  exceeds  the  initial  value  tq.  Once  the  initial  damage  threshold  is  exceeded,  is  updated 
from  Equation  2.2,  dn+i  =  (^(^n+i),  and  a  new  damage  threshold  is  updated  as  follows:  r„+i  =  f„+i. 
Thus  the  damage  threshold  increases  as  damage  accumulates. 

One  functional  form  for  G  suggested  by  Mazaxs  in  Reference  13  is: 


(2.5) 


A  and  B  are  material  constants  which  axe  obtained  from  fits  to  laboratory  test  data.  This  functional 
form  for  G  has  been  implemented  into  the  APTEK  concrete  model.  A  plot  of  G  versus  f  is  shown 
in  Figure  2.1  for  tq  =  14.76  psi”^/^,  A  =  0.8  and  B  =  0.015  psi“^/^.  Note  that  G  =  0  when  f  =  ro 
and  G  approaches  one  as  f  increases. 


2.3  ALTERNATE  FORMULATION. 

The  damage  degradation  model  previously  discussed  in  Section  2.1  models  isotropic  degradation  of 
the  moduli.  This  means  that  both  the  bulk  and  shear  moduli  are  degraded  simultaneously.  However, 
the  preload  data  previously  shown  in  Figure  1.4  suggests  that  bulk  modulus  degradation  may  not 
be  appropriate  for  modeling  the  response  of  concrete  at  low  or  moderate  confining  pressures.  The 
data  suggests  that  daumage  accumulates  with  increasing  deviatoric  strain  or  stress  (uniaxial  strain 
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Figure  2.1.  Increase  in  the  damage  function,  (?,  with  equivalent  strain,  t. 

preload)  but  not  with  increasing  volumetric  strain  or  pressure  (isotropic  compression  preload).  The 
data  indicate  a  substantial  reduction  in  loading  modulus  (Youngs’s  modulus)  following  preloading  in 
uniaxial  strain  to  a  pressure  of  138  MPa.  However,  the  stress-strain  responses  from  the  unconfined 
compression  test,  measured  with  and  without  preloading  in  isotropic  compression,  are  in  reasonable 
agreement.  The  data  do  not  indicate  a  reduction  in  loading  modulus  following  preloading  in  isotropic 
compression  to  a  pressure  of  138  MPa. 

To  account  for  this  behavior,  APTEK  modified^  the  concrete  model  to  eliminate  damage  in  isotropic 
compression,  base  damage  on  the  accumulation  of  distortional  strain  energy,  and  degrad  the  shear 
modulus  only.  Let  Sij  be  the  deviatoric  stress  tensor  defined  as  follows: 

Sij  =  (Tij  -  P8ij  (2.6) 

were  (Tij  is  the  stress  tensor,  P  is  the  pressure,  and  Sij  is  the  Kronecker  delta.  Let  Sij  be  the  effective 
deviatoric  stress  tensor  associated  with  the  undamaged  state.  The  damage  parameter  d  transforms 

^This  alternative  formulation  emd  the  original  isotropic  formulation  axe  separately  maintained. 
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the  effective  deviatoric  stress  tensor  into  the  deviatoric  stress  tensor,  as  follows: 

5o-  =  (1  -  d)Sii  (2-7) 

It  is  assumed  that  the  damage  parameter  d  depends  on  the  history  of  the  deviatoric  strains.  The 
damage  evolution  function  and  the  equivalent  strain,  fd,  are  defined  as  follows: 

d  =  ^(f,)  (2.8) 

fd  =  ^/^  (2.9) 

Here  is  the  undamaged  elastic  distortional  strain  energy,  defined  in  terms  of  the  undamaged 
elastic  strain  energy,  as  follows:  =  '$°  —  Kel  where  e„  is  the  volumetric  strain.  Damage  initiates 

when  fd  exceeds  a  user  specified  threshold,  tq,  as  follows:  fd  >  vq.  Once  the  initial  damage  threshold 
is  exceeded,  d  is  updated  from  Equation  2.8  and  the  new  damage  threshold  is  set  equal  to  the  old 
(previously  exceeded)  value  of  fj. 

2.4  APPLICATION  TO  A  CAP  MODEL. 

We  applied  Simo  and  Ju’s  damage  degradation  formulation  to  a  three-invariant,  viscoplastic  cap 
model.  We  obtained  the  original  two-invariant,  rate-independent  formulation  of  the  cap  model 
from  Pelessone  in  Reference  2  and  implemented  it  into  DYNA3D.  This  formulation  is  reviewed  in 
Section  2.4.1.  We  extended  this  formulation  to  include  the  effect  of  the  third  deviatoric  stress 
invariant  on  geologic  material  strength.  The  three-invariant  model  fits  a  wider  range  of  laboratory 
test  data  than  the  two-invariant  model.  The  three-invariant  formulation  was  developed  by  Rubin 
in  Reference  7  and  is  discussed  in  Section  2.4.2.  Verification  of  the  three-invariant  formulation  is 
given  in  Appendix  C.  We  also  added  rate-dependence  to  the  model  through  implementation  of 
viscoplastic  and  viscodamage  formulations.  These  implementations  allow  us  to  model  an  increase 
in  strength  as  a  function  of  strain  rate.  They  were  formulated  by  Simo  and  Ju  in  Reference  5  and 
are  discussed  in  Section  2.4.4. 

We  also  extended  the  cap  model  to  include  explicit  dependence  of  the  pressure-volumetric  strain 
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(a)  Smooth-cap. 


(b)  Non-smooth  cap. 


Figure  2.3.  Comparison  of  the  smooth-cap  yield  surface  used  in  the  APTEK  concrete  model  with 
a  non-smooth  cap  yield  surface. 


relation  on  the  current  porosity.  The  user  has  the  option  of  requesting  the  explicit  pore-compaction 
model  or  the  original  constant  bulk  modulus  model.  The  explicit  pore  compaction  model  was 
developed  for  soils  and  is  discussed  in  Appendix  B.  The  constant  bulk  modulus  model  was  used  to 
model  concrete  in  all  our  benchmark  and  internal  detonation  calculations. 


2.4.1  Smooth-Cap  Formulation. 

The  cap  model  proposed  by  Pelessone  is  a  two-invariant  plasticity  model  with  a  yield  surface 
defined  by  a  shear  failure  envelope,  a  moveable  hardening  cap,  and  a  tensile  cutoff.  There  is  a 
continuous  and  smooth  (no  corners)  intersection  between  the  failure  envelope  and  hardening  cap, 
which  eliminates  the  numerical  complexity  of  treating  a  compressive  ‘corner’  region  between  the 
failure  surface  and  cap^.  Hence  we  refer  to  Pelessone’s  cap  model  as  a  smooth-cap  model.  Figure  2.2 
shows  a  three-dimensional  view  of  a  two-invariant  cap  model  in  principal  stress  space.  Figure  2.3 
shows  a  two-dimensional  view  in  stress  invariant  space,  along  with  a  non-smooth  yield  surface  for 
comparison.  The  nomenclature  used  in  these  figures  will  be  explained  in  subsequent  paragraphs. 


The  two-invariant  model  is  described  in  terms  of  two  stress  invariants,  Ji  and  Jj,  were  Ji  is  the 
first  invariarit  of  the  stress  tensor: 


Ji  =  ffkk  =  3P  (2.10) 


^One  goal  of  the  cap  model  proposed  by  Pelessone  was  to  improve  the  cdgorithmic  implementation  of  the  Simo- 
Ju-Taylor  cap  model  previously  implemented  in  the  LLNL  version  of  DYNaSd.  In  the  Simo-Ju-Taylor  model,  the  cap 
has  a  horizontal  tangent  that  intersects  the  sloped  failure  surface  in  a  non-smooth  fashion,  i.e.  the  cap  and  feiilure 
surface  have  different  slopes.  Hence  special  treatment  must  be  given  to  the  corner  region  formed  by  two  lines  normal 
to  the  cap  and  failure  surfaces. 
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Figure  2.2.  Two-invariant  cap  model  in  principal  stress  space, 
and  J'2  is  the  second  invariant  of  the  deviatoric  stress  tensor: 

j;  =  \SiiSii  (2.11) 

Here  bars  are  used  over  the  stresses  and  invariants  to  indicate  effective  quantities  associated  with 
the  undamaged  state  in  Equation  2.1. 

The  yield  function  of  the  two-invariant  model  is  described  in  terms  of  the  shear  failure  surface  and 
hardening  cap  functions,  Ff  and  Fc  respectively,  as  follows: 


=  (2.12) 

Multiplying  the  cap  ellipse  function  by  the  failure  surface  function  allows  the  cap  to  take  on  the 
slope  of  the  failure  surface  function  at  their  intersection,  as  discussed  in  subsequent  paragraphs. 
Thus  the  yield  surface  is  described  by  a  single  smoothly  varying  function,  /,  as  previously  shown  in 
Figure  2.3.  Stress  states  inside  the  failure  envelope,  /  <  0,  axe  elastic  states.  Stress  states  on  the 
envelope,  /  =  0,  indicate  that  the  material  is  yielding.  Stress  states  outside  the  envelope,  /  >  0, 
are  not  allowed. 
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The  shear  failure  surface  is  defined  in  terms  of  the  first  stress  invariant  as; 


Ff{Ji)  =  a-  +  eJi 


(2.13) 


The  material  parameters  a,  A,  (3,  and  6  are  evaluated  by  fitting  Equation  2.13  to  experimental  data 
{^fj^  versus  Ji)  from  TXC  tests,  and  then  adjusting  these  parameters  to  account  for  damage,  as 
discussed  in  Appendix  D. 

The  isotropic  hardening  cap  Fc{Ji,  k)  is  a  two-part  function  that  is  either  unity  (thus  reducing  the 
yield  function  to  Fj{Ji)  for  the  shear  failure  surface)  or  an  elliptically  shaped  function  for  the  cap 
with  size  and  position  (hardening)  parameter  k.  The  cap  surface  is  defined  as: 


°  nx(K)-L{K)r 


(2.14) 


This  function  is  displayed  schematically  in  Figure  2.4.  It  is  unity  for  Ji  less  than  L{k)  and  elliptical 
for  L{k)  <  Ji  <  X(k).  The  hardening  parameter  k  that  controls  the  motion  of  the  cap  surface, 
and  L(k)  and  A'(/c)  define  the  geometry  of  the  cap.  The  function  L(k)  is  the  current  value  of  Ji  at 
the  intersection  of  the  cap  and  failure  surfaces.  The  function  A’(k)  is  the  current  value  of  Ji  at  the 
intersection  of  the  cap  and  the  Ji  axis.  It  is  defined  in  terms  of  the  cap  ellipticity  ratio  R,  where  R 
is  the  ratio  of  its  major  and  minor  axes  (parallel  to  the  Ji  to  axes),  as  follows: 


^(/c)  =  L(k)  RFf(L(K)) 


(2.15) 


where  L(k)  is  defined  by: 


T(«)  = 


K  if  K  >  Ko 

Ko  otherwise 


(2.16) 


Here  kq  is  the  value  of  Ji  at  the  initial  intersection  of  the  cap  and  failure  surfaces.  Note  that 
Equation  2.16  restrains  the  cap  from  retracting  past  its  initial  location  at  k).  This  restraint  is 
arbitrary  and  the  implementation  can  be  readily  modified  by  the  user  to  allow  the  cap  to  retract 
past  its  initial  location,  should  such  a  modification  be  warranted  by  data.  In  fact,  the  user  has  the 
option  of  not  allowing  the  cap  to  contract  at  all,  which  is  often  the  desired  option  for  modeling 
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Figure  2.4.  Two-paxt  function  used  for  the  cap. 


rock. 

When  plastic  volume  compaction  or  expansion  occurs,  the  evolution  of  the  cap’s  motion  is  defined 
by  the  hardening  rule: 

ej  =  W  (l  —  (2-17) 

where  ej  is  the  plastic  volumetric  strain,  W  is  the  maximum  plastic  volumetric  strain,  Xo  =  X(ko) 
is  the  initial  abscissa  intercept  of  the  cap  surface,  and  J?i  and  D2  are  material  parameters.  The 
four  parameters  of  the  hardening  rule  are  obtained  from  the  pressure- volumetric  strain  curve  from 
isotropic  compression  tests.  The  cap  aspect  ratio  R  is  typically  evaluated  from  a  fit  to  the  stress 
path  from  uniaxial  strain  data.  A  detailed  description  of  how  to  fit  the  model  to  concrete  data  is 
given  in  Appendix  D  for  the  benchmark  series  of  calculations. 

2.4.2  Three  Invariant  Extension. 

We  extended  the  two-invariant  smooth-cap  model  to  include  the  well  known  dependence  of  geologic 
materi2d  behavior  on  the  third  invariant  of  the  deviatoric  stresses,  J3.  This  extension  uses  a  formu¬ 
lation  published  by  Prof.  Rubin  in  Reference  7.  The  three-invariant  model  describes  more  complete 
geologic  material  behavior  from  a  wider  range  of  laboratory  test  data  than  the  two-invanant  model. 

Standard  laboratory  tests  for  measuring  failure  curves  of  geologic  specimens  include  triaxial  com¬ 
pression  (tXC)  and  triaxial  extension  (tXE)  tests.  Such  tests  are  conducted  on  cylindrical  specimens 
and  begin  with  hydrostatic  compression,  t.e.  the  axial  compressive  stress,  ax,  is  equal  to  the  radial 
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Figure  2.5.  Schematic  failure  curves  for  triaxial  compression  (txc),  a  deviatoric  state  of  torsion 
(tor),  and  triaxial  extension  (txe). 

compressive  stress,  <7,.  For  TXC  tests,  the  magnitude  of  the  compressive  stress  is  then  quasi-statically 
increased  until  the  specimen  fails.  For  TXE  tests,  the  magnitude  of  the  compressive  stress  is  quasi- 
staticeilly  decreased  until  the  specimen  fails.  Typical  failure  curves,  shown  in  Figure  2.5,  are  plotted 
as  a  stress  difference  quantity,  =  {cTx  —  <Tr)/\/3,  versus  pressure,  Ji/3.  It  is  well  known  that 
geologic  materials  fail  at  lower  values  of  for  TXE  tests  than  for  TXC  tests  conducted  at  the 
same  pressure.  The  failure  surface  functions  Fj,  Qi,  and  shown  in  this  figure  are  discussed  in 
subsequent  paragraphs. 

A  two-invariant  model  cannot  be  simultaneously  fit  to  measured  failure  curves  from  both  TXC  and 
TXE  tests.  Two-invariant  models  are  typically  fit  to  TXC  data,  although  the  choice  of  which  data 
(txc  or  txe)  to  fit  depends  on  the  analyst’s  judgement  and  the  availability  of  data.  The  advantage 
of  the  three-invariant  model  proposed  by  Rubin  is  that  it  can  be  simultaneously  fit  to  measured 
TXC  and  TXE  failure  curves,  as  well  as  to  a  failure  curve  for  a  deviatoric  state  of  torsion  (tor). 

Figure  2.6  shows  a  three-dimensional  view  of  a  three-invariant  cap  model  in  principal  stress  space. 
Although  a  straight-line  Mohr-Coulomb  failure  envelope  is  shown  in  Figure  2.6,  our  model  is  more 
general  than  shown  because  it  adlows  the  user  to  specify  a  nonlinear  failure  curve  to  describe  the 
pressure  dependence  of  the  TXC  curve  (refer  back  to  Equation  2.13).  The  difference  between  the 
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Figure  2.6.  Three-invajictnt  cap  model  in  principal  stress  space. 

two-invariant  cap  model  previously  shown  in  Figure  2.2  and  the  three-invariant  cap  model  shown  in 
Figure  2.6  is  the  shape  of  the  yield  function  in  the  octahedral  plane^.  For  the  two-invariant  model, 
the  yield  function  forms  a  circle  in  the  octahedral  plane.  For  the  three-invaxiant  model,  the  yield 
function  forms  a  hexagon  (regular,  semi-regular  as  shown,  or  irregular)  in  the  octahedral  plane. 

Figure  2.7  depicts  three  common  failure  surfaces  in  the  octahedral  plane.  The  two-invariant  model 
corresponds  to  the  the  circular  surfaces,  e.g.  Mises-Schleicher  or  Drucker-Prager.  The  three- 
invariant  model  corresponds  to  the  Mohr-Coulomb  semi-regulax  hexagon.  It  is  obtained  fits  to 
both  TXC  (point  A)  and  TXE  (point  C)  data.  Although  a  straight-line  is  fit  between  the  TXC  and 
TXE  states,  the  Rubin  three-invariant  formulation  is  general  enough  that  the  yield  function  may 
form  an  irregular  hexagon-like  shape  in  which  each  of  the  six  sides  is  quadratic  (rather  than  linear 
as  shown)  between  the  TXC  and  TXE  states.  Such  a  shape  is  obtained  from  fits  to  TXC,  TOR,  and 
TXE  data. 

The  two-invariant  model  was  extended  to  three  invariants  by  modifying  the  yield  function  in  Equa- 

^The  octahedral  plane  is  a  plane  whose  normal  is  the  hydrostatic  a.xis  and  makes  equal  angles  with  each  of  the 
principal  axes  of  stress. 
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Figure  2.7.  Three  common  failure  surfaces  in  the  octahedral  plane. 


tion  2.12,  as  follows: 

/( J„  j;,  «)  =  /,-  Tl^F]F,  (2.18) 

where  Ji)  is  a  scaling  function  proposed  by  Rubin  in  Reference  7.  The  angle  ^  is  confined  to 
the  range  — 7r/6</3<7r/6,  and  is  related  to  the  invariants  ^xid  Jg  as  follows: 

sin3^  =  j3  =  2^^  (2.19) 

in  which  Jg  is  the  determinant  of  the  deviatoric  stress  tensor: 

j'z  =  (2.20) 

A  A 

Jz  is  a  normalized  invariant  which  remains  in  the  range  —  1  <  J3  <  1.  For  the  standard  laboratory 
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tests  just  discussed,  the  values^  of  ^  and  J3  axe: 


=  I,  J3  =  1  for  TXC 

0 

/9  rr  0,  J3  =  0  for  TOR  (2.21) 

/3  =  -^,  J3  =  -l  forTXE 
0 

The  scaling  function  K  in  Equation  2.18  determines  the  strength  of  the  material  for  any  state 
of  stress  relative  to  the  strength  for  TXC.  This  is  demonstrated  by  referring  back  to  Figure  2.5. 
Three  functions  of  pressure  axe  indicated,  denoted  Fj,  Qi,  and  Q2-  The  function  Ff  is  fit  to  the 
failure  curve  measured  from  TXC  tests.  The  product  functions  QiFj  and  Q^Ff  are  fit  to  failure 
curves  measured  from  TOR  and  TXE  tests.  Rubin  develops  an  analytical  expression  for  TZ  in  terms  of 
experimental  values  for  Q\  and  Q2  as  functions  of  pressure.  This  expression  in  given  in  Equation  A.2 
of  Appendix  A.  The  scaling  function  TZ  takes  on  positive  values  less  than  or  equal  to  one,  and  thus 
it  scales  the  failure  curve  Ff  for  stress  states  other  than  TXC. 

We  choose  to  implement  the  following  functional  form  for  the  scaling  fimctions  Qi  and  Q2- 

Qi  =  ai  -  Ji  (2.22) 

Q2  =  012  —  ^26  +  O2J1  (2.23) 

The  eight  material  parameters  (ai,Ai,^i,^i,  and  012^21^21^2)  evaluated  by  fitting  the  product 
functions  QiFf  2ind  Q2Ff  to  experimental  data  {\fj2  versus  Ji)  from  TOR  and  TXE  tests,  and  then 
adjusting  the  peirameters  to  account  for  damage.  Minimum  values  axe  Qi  =  0.57735  and  Q2  =  0.5. 
These  minimum  values  define  a  triangle  in  the  octahedral  plane.  Maximum  values  are  Qi  =  1.0  and 
Q2  =  1.0.  These  maximum  values  define  a  circle  in  the  octahedral  plane.  The  Rubin  three-invariant 
formulation  allows  the  shape  of  the  yield  function  in  the  deviatoric  plane  to  transition  with  pressure 
from  triangular,  to  irregular  hexagonal,  to  circtilax. 

The  yield  function  provided  by  Equation  2.18  (with  Fc  =  1)  reduces  to  the  classical  criteria  of 
^Stresses  are  assumed  positive  in  compression. 
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Mohr-Coulomb  aad  maximum-octaiiedral-sheax  stress  for  specific  values  of  Q\  and  Q2,  and  specific 
functional  forms  for  Fj.  The  reduced  forms  of  %  axe  given  in  Appendix  A.  Verification  of  the  three- 
invariant  model  for  constant  (pressure-independent)  values  for  Qi  and  Q2  is  given  in  Appendix  C. 
Reduction  of  the  Rubin  model  to  a  Wiliam- Warnke  model  is  also  given  in  Appendix  A.  The  Willam- 
Warnke  model  (Reference  15)  is  a  commonly  used  model  in  which  the  shape  of  the  yield  surface  in 
the  deviatoric  plane  transitions  with  pressure:  the  Wiliam- Waxnke  model  is  a  subset  of  the  Rubin 
model. 


2.4.3  Implementation  Aspects. 

The  implementation  of  the  damaged  degradation  formulation  into  finite  elements  models,  such  as 
the  three-invaxiant  cap  model,  is  particularly  simple  because  it  makes  use  of  an  additive  split  of 
the  the  stress  tensor  into  elastic-damage  and  plastic  corrector  parts.  Damage  evolution  and  plastic 
flow  axe  treated  as  separate  processes.  The  algorithmic  formulation  proceeds  in  two  phases:  an 
elastic-damage  prediction  in  which  the  parameter  d  is  updated,  and  a  plastic  return  mapping  phase 
in  which  the  effective  stresses,  a-,j,  are  updated. 

Damage  Parameter  Update.  Let  n  denote  the  nth  time  step  in  the  finite  element  analysis.  The 
first  step  in  updating  the  damage  parameter  from  to  dn+i  is  to  update  the  total  strain  tensor, 
e„,  from  the  incremental  strain  tensor,  Ac,  as  follows:  tn^i  =  Cn  +  Ae.  Then  one  computes  the 
equivalent  strain  according  Equation  2.3  or  Equation  2.9  as  follows: 

f„+i  =  y^2$o(e„+i)  (isotropic)  (2.24) 

The  updated  damage  parameter  and  damage  threshold  r„+i  axe  then  given  by: 


1  dn  if  Tn+i  -  r„  <  0 

dn+1  =  < 

1  G(f„+i)  otherwise 

(2.25) 

r„+i  =  max(r„,f„+i) 

(2.26) 
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Effective  Stress  Update.  The  effective  stress  components,  dij  in  Equation  2.1,  are  those  stress 
components  obtained  from  the  elasto-plastic  cap  model  without  consideration  of  damage.  The 
recovery  of  the  effective  stress  components  at  the  beginning  of  the  time  step  is  straight  forward, 
t.e.  o-„  =  o-n/(l  —  dn).  Plastic  response  is  formulated  in  effective  stress  space  holding  d„  constant. 

The  effective  stress  update  to  determine  5-n+i  requires  several  steps.  The  first  step  is  to  check  for 
yielding.  This  is  done  by  by  temporarily  updating  the  effective  stress  components  from  the  incre¬ 
mental  strain  rates  by  assuming  that  the  entire  strain  increment  is  elastic.  These  updated  stresses 
are  called  the  elastic  trial  stresses,  denoted  .  The  value  of  the  yield  function  in  Equation  2.18 
is  valuated  from  the  elastic  trial  stresses.  The  new  stress  state  is  an  elastic  stress  state  if  /  <  0  and 
a  plastic  stress  state  is  /  >  0. 


If  the  stress  state  lies  inside  the  failure  surface  (/  <  0)  then  the  elastic  stress  update  is  trivial, 
t.e.  If  the  stress  state  lies  outside  the  failure  surface  (/  >  0),  then  a  return  mapping 

algorithm,  or  plasticity  algorithm,  returns  the  stress  state  back  to  the  failure  surface  so  that  /  =  0. 
This  is  done  by  enforcement  of  the  so-called  plastic  consistency  condition,  /  =  0,  which  may  be 
expressed  as: 


/  =  M-^i:  + 


da, 


*3 


I'- 


Here  the  superscript  dot  indicates  differentiation  with  respect  to  time. 


(2.27) 


Enforcement  of  the  consistency  condition  requires  several  assumptions.  The  first  of  these  assump¬ 
tions  is  that  the  total  strain  rate  tensor,  c,j,  cam  be  partitioned  into  elastic  and  plastic  components: 


^3  = 


(2.28) 


The  goal  of  the  plasticity  algorithm  is  to  determine  the  plastic  strain  rate,  thus  identifying  the 
partitioning. 

The  second  zissumption  is  that  the  plastic  strain  increment  is  normal  to  the  yield  surface: 


(2.29) 
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where  A  is  a  proportionality  constant  known  as  the  consistency  parameter.  This  assumption  is 
known  as  an  associated  flow  rule,  or  normality  condition.  Use  of  a  potential  function  other  than 
the  yield  function  in  Equation  2.29  results  in  a  nonassociated  flow  rule.  Recent  studies  reported 
in  Reference  16  suggest  that  rate-independent  models  with  nonassociated  flow  lead  to  spurious 
(non-unique)  dynamic  solutions,  so  only  associated  flow  is  implemented  in  the  present  model. 

The  third  assumption  defines  the  evolution  relation  for  the  moveable  cap,  known  as  the  hardening 
rule: 

k  =  \h{dij,K)  (2.30) 

where  h  is  the  plastic  hardening  law.  For  the  smooth-cap  model,  we  readily  define  h  by  equating 
the  plastic  volumetric  strain  rate,  ej,  derived  from  the  flow  rule  in  Equation  2.29  to  that  obtained 
from  a  first  order  Taylor  series  expansion  of  the  hardening  rule  in  Equation  2.17: 

M^.,«)  =  3(|i)/(f)  (2.31) 

One  additional  assumption  is  the  form  of  the  constitutive  relation.  For  concrete,  we  assume  Hooke’s 

Law  for  an  isotropic  material: 

-  cfj)  (2.32) 

where  the  elastic  constitutive  tensor  Cijki  from  Reference  17  is  given  by: 

Ciiki  =  {K-  ‘^G)6ii6ki  +  G{8ikSii  +  SuS^k)  (2.33) 

For  soil,  another  form  of  the  constitutive  relation  is  the  explicit  pore  compaction  model  discussed 
in  Appendix  B. 

The  solution  of  the  consistency  condition  in  Equation  2.27  determines  A,  which  in  turn  determines 
the  partitioning  of  the  total  strjiin  rate  into  elastic  cind  plastic  components.  The  A  solution  is 
obtained  from  Equation  2.27  through  substitution  of  the  yield  function  in  Equation  2.18,  the  flow 
rule  in  Equation  2.29,  the  hardening  rule  in  Equation  2.31,  and  the  assumed  constitutive  relation 
in  Equation  2.32. 
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Details  of  the  solution  procedure  are  given  in  Reference  18  for  the  special  case  of  'R,  independent  of 
Ji  and  a  linear  elasticity  tensor,  i.t.  Hookes  law  with  constant  undamaged  moduli  K  and  G.  For 
this  special  case,  the  expression  for  A  is: 


(2.34) 


with  the  compact  notation 

=  ^  (2.35) 

Details  of  the  solution  procedure  for  more  sophisticated  constitutive  models,  such  as  the  explicit 
pore  compaction  model  discussed  in  Appendix  B,  have  not  been  summarized  due  to  a  shift  in 
interest  on  the  CWE  program  from  modeling  soils  to  modeling  concrete.  The  effective  bulk  modulus 
used  in  the  explicit  pore  compaction  model  varies  with  k.  Such  dependence  of  the  modulus  on  the 
hardening  parameter  must  be  taken  into  account  in  the  solution  procedure,  but  details  are  beyond 
the  scope  of  this  report. 


2.4.4  Rate-Dependent  Extensions. 

Some  existing  experimental  data  shows  that  the  peak:  stress  attained  during  direct  pull  and  uncon¬ 
fined  compression  tests  is  sensitive  to  the  rate  of  loading  or  the  strain  rate.  An  example  of  such  data 
from  Reference  19  is  reproduced  in  Figure  2.8.  To  accommodate  such  strain-rate  sensitivity,  APTEK 
implemented  a  rate- dependent  model  which  requires  one  additional  parameter,  denoted  by  r\.  This 
parameter  is  called  the  fluidity  coefficient,  and  the  rate-dependent  model  is  called  a  viscoplastic 
model.  The  original  constitutive  formulation  is  attributed  to  Duvaut  and  Lions  in  Reference  14. 

Simo  and  Ju  in  Reference  4  suggest  that  the  amount  of  microcracking  (damage)  exhibited  by  some 
materials  at  a  psui^icular  strain  level  is  retarded  at  higher  strain  rates,  thus  decreasing  the  non¬ 
linearity  of  the  stress-strain  curve  at  higher  strain  rates.  To  model  such  behavior,  they  suggest 
implementation  of  a  rate-sensitive  viscodamage  model  which  requires  one  additional  damage  fluid¬ 
ity  coefficient,  denoted  by  r^.  The  viscoplastic  and  viscodamage  extensions  are  discussed  in  the 
following  paragraphs. 
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Figure  2.8.  Measured  rate-dependence  of  the  peak  stress  attained  in  direct  pull  and  unconfined 
compression  tests. 

Viscoplastic  Extension.  Simo  et.al.  in  Reference  12  postulate  a  three-dimensional  generalization 
of  the  Duvaut-Lions  viscoplastic  strain  rate  formulation  as  follows: 


V 


(2.36) 


where  is  the  viscoplastic  strain  rate  tensor,  C  is  the  elasticity  tensor,  and  a  and  5"  are  the  viscid 
and  inviscid  stress  tensors,  respectively.  For  simplicity,  we  will  use  bold  symbols  in  this  section  to 
indicate  tensors,  i.e.  er  =  ay,  and  the  symbol  to  indicate  tensor  contraction. 


Simo  et.al.  extended  the  Duvaut-Lions  formtilation  to  multi-surface  cap  models  by  developing  a 
closed-form,  first  order  aiccurate  update  solution  which  accommodates  an  arbitrary  number  of  in¬ 
ternal  state  vairiables,  such  as  the  hardening  vau-iable  k.  The  update  solution  requires  the  corre¬ 
sponding  solution  of  the  inviscid  elasto-plastic  problem,  which  we  denote  as  d„+i  and  k„+i.  The 
viscoplastic  update  solution  also  requires  that  the  elastic  trial  stress,  be  saved. 

The  update  algorithm  for  the  stress  tensor  is  obtained  by  direct  application  of  an  implicit  backward 
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Simo  et.al.  also  derives  a  general  update  algorithm  for  the  internal  variables, 
that  algorithm  for  the  hardening  variable  k  is: 

/c„  +  At/ri 

l  +  Ai/, 

Note  that  the  inviscid  and  elastic  solutions  can  be  obtained  with  these  update  algorithms.  The 
inviscid  solution  (o-„+i  =  5-„+i)  is  obtained  from  Equation  2.42  as  t]  0.  This  situation  corresponds 
to  the  rate-independent  solution.  The  elastic  solution  (o-n+i  =  (t^^)  is  obtained  as  7  — »  oo.  This 
situation  corresponds  to  the  absence  of  plastic  flow.  At  each  timestep,  the  viscoplastic  solution  is 
bounded  between  the  current  rate-independent  elasto-plastic  solution  and  the  elastic  trial  solution. 

Numerous  single  element  simulations  were  preformed  for  unconfined  compression  to  determine  the 
numerical  relationship  between  strength  enhancement,  strMn  rate,  and  fluidity  parameter.  By 
strength  enhancement  we  mean  the  ratio  of  the  peak  stress  attained  with  viscoplasticity  {t]  >  0) 
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relative  to  the  peak  stress  attained  without  viscoplasticity  (7/  =  0).  We  will  refer  to  this  ratio  as  the 
strength  enhancement  factor,  SEF.  Our  simulations  indicate  that  the  calculated  strength  enhance¬ 
ment  factor  is  nonlinearly  related  to  the  product  of  the  strain  rate  times  the  fluidity  parameter,  erj, 
as  shown  by  the  diamond  symbols  in  Figure  2.9  for  our  WSMR-5  concrete  model.  The  nonlinear 
relationship  is  material  dependent  and  time-step  independent.  Note  that  Figure  2.9  provides  a 
calculated  SEF  factor,  whereas  Figure  2.8  provides  a  measured  SEF  factor. 


If  a  calculation  is  to  be  run  at  a  constant  strain  rate,  then  the  user  can  readily  select  a  fluidity 
parameter  value  match  the  measured  strength  enhancement  factor  at  that  particular  strain  rate. 
However,  the  current  strain  rate  in  most  calculations  varies  on  an  element-by-element  and  time- 
step  by  time-step  basis.  For  such  a  case,  we  propose  a  method  of  internally  calculating  the  fluidity 
parameter  based  on  a  user  supplied  function  for  the  measured  strength  enhancement  factor. 


Suppose  the  measured  strength  enhancement  vs.  strain  rate  is  formulated  as  follows: 


SEF{e) 


1.29  —  .059  *  log(c)  for  e  >  1.0e~® 
0  otherwise 


(2.44) 


where  e  is  the  effective  strain  rate: 


and  is  the  deviatoric  strain: 


(2.45) 


(2.46) 


This  function  is  plotted  in  Figure  2.10  and  can  be  obtained  from  a  fit  to  rate-dependent  data, 
such  as  the  data  previously  shown  in  Figure  2.8.  Equation  2.44  is  the  benchmark  specification  for 
WSMR-5  concrete. 


Internally  calculating  the  fluidity  parameter  is  a  two-step  process.  The  first  step  is  to  calculate  the 
specified  SEF  from  Equation  2.44  (or  Figure  2.10)  given  i  in  the  element  of  interest.  The  second 
step  is  to  back-calculate  t)  from  the  specfied  SEF  and  e.  This  is  done  by  fitting  a  functional  form  to 
the  calculated  SEF  versus  et)  values  previously  shown  in  Figure  2.9.  Then  the  product  er]  is  back- 
calculated  from  this  functional  form,  or  equivcilently,  estimated  from  the  graphical  representation 


26 


Figure  2.9.  Calculated  rate-dependent  strength  enhancement  factor  established  from  unconfined 
compression  simulations  for  WSMR-5  concrete. 


Figure  2.10.  Benchmark  specification  for  strength  enhancement  vs.  strain  rate  estimated  from  mea¬ 
sured  data. 
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in  Figure  2.9.  The  fluidity  pajameter  t)  is  readily  derived  from  the  product  et]  because  the  strain 
rate  e  is  known. 


Viscodamage  Extension.  Simo  and  Ju  also  developed  a  closed-form  update  algorithm  in  Ref¬ 
erence  5  for  linear  viscous  damage.  Suppose  f„+i  >  r„  so  that  damage  accumulation  is  predicted. 
The  update  algorithms  for  the  damage  parameter  and  damage  threshold  are: 


dn+i 


=  d„  d-  Ar 


dGjfn) 

df 


^n+i  =  r„  -f  Ar 


(2.47) 

(2.48) 


where  Ar  =  A/i(T„+i  —  r„). 
as  follows: 


The  parameter  A/i  is  derived  from  the  damage  fluidity  coefficient,  rjd, 


Afi  = 


Atlr)d 

1  -f  At/rjd 


(2.49) 


Both  the  instantaneous  elasticity  and  inviscid  damage  solutions  are  obtained  as  limiting  cases  of  the 
rate-dependent  formulation.  The  instantaneous  elasticity  solution  (r„+i  =  r„  and  dn+i  =  =  0) 

is  obtained  as  rfd  oo.  The  inviscid  solution  (r„+i  =  fn+i  and  =  4  -h  Af^)  is  obtained 
from  Equation  2.48  as  t/j  0.  Hence  the  expansion  of  the  damage  surface  (r„  <  r„+i  <  f„+i)  is 
properly  bounded  between  the  instantaneous  elasticity  and  the  inviscid  damage  limit. 


2.4.5  Single  Element  Elasto-Plcistic  Damage  Calculation. 

Results  of  a  uniaxial  compressive  stress  calculation  under  cyclic  loading,  which  includes  damage  and 
plasticity,  are  given  in  Figure  2.11.  The  calculated  stress-strain  response  captures  the  essential 
features  of  the  measured  response  previously  shown  in  Figure  1.3.  No  attempt  was  made  to  fit  our 
model  to  this  data  because  not  enough  information  was  available. 

One  feature  of  the  damage  model  demonstrated  in  Figure  2.11  is  strain-softening,  which  is  a  decrease 
in  strength  during  progressive  straining  after  a  peak  strength  value  is  reached.  Another  feature  is 
stiffness  degradation.  The  average  slopes  of  the  loading/unloading  curves  decrease  with  increased 
straining,  which  indicates  progressive  degradation  of  the  elastic  moduli.  The  damage  paxcimeter  d 
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Figure  2.11.  Uniaxial  compressive  stress  calculation  showing  softening  and  degradation  of  the  elastic 
moduli  with  increased  straining. 

remains  constant  during  unloading  and  reloading  because  unloading  reduces  the  equivalent  strain 
below  the  damage  threshold.  Also  note  that  permanent  plastic  strains  occur  upon  unloading. 

This  example  calculation  is  intended  to  demonstrate  pertinent  features  of  the  damage  model  but 
the  model  parameters  were  not  obtained  from  a  fit  to  laboratory  test  data.  Fits  of  the  model  to 
WSMR-5  concrete  data  axe  given  in  Section  7  for  benchmark  calculations. 
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SECTION  3 


MODELING  DAMAGE  IN  THE  BRITTLE  REGIME 


Of  great  importance  to  the  OWE  program  is  the  behavior  of  concrete  in  the  low  confining  pressure 
and  tensile  (brittle)  regimes.  We  define  the  brittle  regime  as  any  state  of  stress  in  which  the 
pressure  is  tensile.  Pelesonne’s  original  cap  model  included  implementation  of  a  tensile  pressure 
cutoff  model.  This  model  set  the  stress  equal  to  zero  (over  one  time  step)  once  the  pressure  cutoff 
value  was  exceeded.  We  used  this  model  for  the  Round  A  and  B  benchmarks,  and  for  analysis  of 
a  Dipole  East  structure  (deint-3).  However,  our  DEINT-3  calculation  did  not  correlate  well  with 
the  test  data  because  the  one-step  reduction  in  tensile  strength  was  too  severe.  We  obtained  better 
correlations  in  subsequent  structural  response  calculations  (Large  Test  Structure  1  and  the  Precision 
Wall  Test)  by  controling  the  rate  of  softening  in  the  brittle  regime. 

To  improve  the  model  in  this  regime,  and  to  meet  the  benchmark  specifications,  we  eliminated  the 
tensile  cutoff  portion  of  the  model  and  implemented  an  elasto-damage  model  for  brittle  response 
(plasticity  is  neglected).  Our  damage  model  in  the  brittle  regime  (P  <  0)  was  specifically  formu¬ 
lated  to  meet  all  of  RDa’s  peak  stress  specifications  for  the  Round  D  and  E  benchmark  series  of 
calculations.  These  specifications  are  given  in  Table  3.1  and  indicate  that  the  maximum  principal 
stress  cannot  exceed  a  peak  stress  of  3.2  MPa.  WSMR-5  test  data  in  the  brittle  regime  was  limited  to 
the  axial  stress-strain  response  for  direct  pull,  measured  with  and  without  the  effects  of  preloading 
in  uniaxial  strain  and  isotropic  compression.  This  data  was  previously  shown  in  Figure  1.4.  Hence 
the  peak  stress  specifications  listed  in  Table  3.1,  other  than  direct  pull,  were  subjectively  selected 
by  RDA  to  provide  a  complete  (three-invariant)  failure  surface. 

Table  3.1.  RDA’s  benchmark  specifications  in  the  brittle  regime. 


Simulation 

Peak  Stress  (MPa) 

Direct  PuU 

3.2 

Pure  Shear 

3.2 

Equal  Biaxial  Tension 

3.2 

Equal  Triaxial  Tension 

3.2 
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3.1  STRESS-BASED  FORMULATION. 


Simo  and  Ju  formulate  an  anisotropic  model  for  brittle  response  in  Reference  4  which  requires 
second-  and  fourth-order  tensor  damage  variables^.  However,  we  chose  to  loosely  base  our  brittle 
damage  model  on  their  isotropic  formulation  in  Equation  2.1,  due  to  its  efficiency  and  simplicity 
through  use  of  a  scalar  damage  variable.  In  addition,  we  choose  to  neglect  plastic  response  in 
the  brittle  regime,  due  to  the  limited  availability  of  stress-strain  data.  In  particular,  no  data  was 
provided  on  lateral  strains  or  post-peah  axial  softening  response.  Neglecting  plasticity  makes  the 
model  more  run-time  efficient,  which  is  an  advantage  when  performing  complex  analyses,  such  the 
the  Dipole  East  series  of  calculations.  At  the  same  time,  we  can  readily  fit  the  elastic-damage  model 
to  the  limited  set  of  WSMR-5  data  and  meet  all  of  RDA’s  peak  stress  specifications. 

To  model  the  elastic-damage  response  in  the  brittle  regime,  the  user  specifies  the  peak  stress,  which 
in  this  case  is  Cmax  =  3.2  MPa,  and  two  parameters  {Ab  and  Bb)  which  regulate  the  post-peak 
softening  response  according  to  Equation  2.5.  Damage  initiates  when  an  energy  norm,  n,  exceeds 
iin  initial  threshold,  tq. 

Unlike  the  pressure  independent  threshold  supplied  by  the  user  for  the  ductile  model,  the  threshold 
for  the  brittle  model  is  pressure  dependent.  This  is  necessary  to  meet  the  benchmark  specifications. 
The  initial  damage  threshold  is  calculated  firom  the  user-supplied  peak  stress  specification  and 
undamaged  elastic  moduli,  as  follows: 


where  It  is  the  value  of  the  Rubin  function  for  the  three-invariant  formulation,  which  is  7^  =  0.5  for 
the  benchmark  specifications  in  the  brittle  regime,  and  Yt  is  a  function  of  pressure  which  describes 
the  failure  (peak  stress)  surface,  as  follows: 

YT  =  ^{P  +  a^^)  for  -ar„^<P<0  (3.2) 

^Theii  elastic-damage  moduli  tensor  is  non-symmetric  for  the  anisotropic  model,  and  thus  damage  is  directional 
in  nature.  The  elastic-d2unage  moduli  tensor  is  symmetric  for  the  isotropic  model. 
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The  pressure  dependent  threshold  in  Equation  3.1  is  a  departure  from  the  work  of  Simo  and  Ju.  A 
simpler  approach  would  be  to  specify  a  pressure-independent  threshold,  f.e,  a  constant  value  for  tq, 
as  discussed  in  Section  3.2.  Although  the  WSMR-5  direct  pull  test  data  could  readily  be  modeled 
with  a  pressure-independent  threshold,  we  could  not  meet  the  benchmark  specifications. 

The  current  value  of  the  energy  norm  is  calculated  from  the  effective  stress  invariants,  Jj 
and  the  undamaged  elastic  moduli,  as  follows; 

Once  the  initial  damage  threshold  is  exceeded,  i.e.  fj,  >  tq,  then  the  updated  value  of  the  scaler 
damage  parameter  is  calculated  from  Equation  2.5  (with  A  =  Ab,  B  =  Bb,  and  f  =  fb),  the  elastic- 
damage  stresses  are  updated  from  Equation  2.1,  and  the  new  threshold  for  damage  is  increased  to 
the  current  value  of  fj,.  The  algorithmic  update  is  given  by  Equations  2.25  and  2.26  in  which  f„+i 
is  replaced  by  7^  „+i.  Note  that  the  effective  stresses  in  Equation  2.1  are  simply  the  linear  elastic 
stresses. 

Our  brittle  damage  criterion  {fb  >  tq)  is  a  stress-based  criterion,  in  which  damage  is  linked  to  the 
history  of  the  undamaged  stress  invariants.  We  could  just  as  easily  have  implemented  a  strain-based 
criterion  in  which  damage  accumulates  according  to  an  energy  norm  of  the  strain  tensor.  Because 
we  are  not  modeling  plasticity  in  the  brittle  regime,  the  strain-based  approach  is  identical  to  the 
stress-based  approach  for  the  direct  pull  test  simulations  given  in  Figure  3.1.  However,  the  two 
approaches  are  not  identical  for  cases  in  which  samples  are  preloaded  into  the  ductile  {P  >  0) 
regime,  causing  plastic  response  to  occur  before  loading  in  the  brittle  {P  <  0)  regime.  A  pressure- 
independent  strciin-based  model  was  used  in  the  Round  C  benchmarks  and  is  briefly  discussed  in 
the  next  section.  Preload  effects,  and  the  applicability  of  strain-based  versus  stress-based  models, 
will  be  examined  in  a  future  effort. 
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3.2  STRAIN-BASED  FORMULATION. 


The  brittle  damage  formulation  previously  discussed  in  Section  3.1  is  a  stress-based,  isotropic  formu¬ 
lation.  An  alternate,  but  not  necessarily  equivalent  approach,  is  to  use  a  strain-based  formulation. 
Our  alternate  approach  is  based  on  a  simplification  of  Simo  and  Ju’s  anisotropic  formulation,  in 
which  we  degrad  K  but  not  G.  This  alternate  formulation  was  used  in  the  Round  C  benchmark 
calctilations  prior  to  RDA’s  specification  of  a  maximum  principal  stress  criterion. 

To  model  damage  in  the  brittle  regime,  the  user  specifies  an  initial  damage  threshold,  r©,  and  two 
parameters  (A  =  Aj  and  B  =  Bb)  which  regulate  the  post-peak  softening  response  according  to 
Equation  2.5.  The  brittle  damage  threshold  is  pressure-independent,  and  damage  initiates  when 
an  energy  norm,  n,  exceeds  the  initial  threshold,  tq.  The  energy  norm  is  based  on  the  tensile 
volumetric  strain,  e^,  as  follows:  n  =  Hence  damage  initiates  and  accumulates  when  the 

volumetric  strain  is  tensile,  rather  than  when  the  pressure  is  tensile. 

This  formulation  predicts  damage  accumulation  with  dilation  for  TXC  simulations  with  low  or  no 
confining  pressure.  In  such  a  case,  our  implementation  compares  the  damage  accumulated  by  brittle 
and  ductile  models  and  selects  the  maximum  damage  parameter  as  the  operative  parameter. 

3.3  RATE-DEPENDENT  EXTENSIONS. 

One  method  of  modeling  rate  dependence  of  the  strength  is  the  viscodamage  extension  previously 
described  in  Section  2.4.4.  Another  method  of  modeling  rate  dependence  of  the  strength  is  to 
scale  the  initial  damage  threshold  previously  given  in  Equation  3.1  as  a  function  of  strain  rate. 
This  is  a  simple  two-step  process.  First,  the  strength  enhancement  factor  is  calculated  from  Equa¬ 
tion  2.44  given  the  effective  strain  rate  in  the  particular  element  of  interest.  Second,  the  initial 
damage  threshold,  ro  is  scaled  by  the  strength  enhancement  factor.  The  viscoplastic  extension  is 
not  applicable  to  the  brittle  regime  because  we  neglected  plasticity. 

Single  element  simulations,  with  and  without  rate-dependence  of  the  strength,  are  given  in  Figure  3.1 
for  the  direct  pull  benchmark  specification.  Rate-dependence  was  modeled  by  scaling  the  initial 
damage  threshold.  We  used  initial  undamaged  moduli  oi  K  =■  15520  Mpa  and  G  =  13890  Mpa  for 
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Figure  3.1.  Direct  pull  test  simulations. 

these  simulations,  which  are  linear  to  the  peak  stress.  We  choose  softening  parameters  of  A  =  1.0 
and  B  =  300  MPa^/^  to  provide  rapid  degradations  in  strength.  One  approach  in  selecting  the 
softening  parameters  is  to  relate  the  area  under  the  calculated  stress-displacement  curve  to  the 
measured  fracture  energy  of  concrete,  Gy,  which  is  often  times  available  in  the  literature. 


3.4  MODELING  PRE-PEAK  DEGRADATION. 

In  addition  to  modeling  post-peak  degradation,  the  user  has  the  option  of  modeling  pre-peak  degra¬ 
dation,  or  hardening,  which  is  evident  in  the  WSMR-5  direct  pull  test  data  previously  shown  in 
Figure  1.4a.  Pre-peak  hardening  can  be  modeled  by  appropriate  selection  of  the  maximum  stress 
parameter,  ir^ax,  and  the  softening  parameters,  Ab  and  B^.  A  comparison  of  three  direct  pull  sim¬ 
ulations  is  shown  in  Figure  3.2.  The  solid  line  simulations  use  a  peak  stress  of  <7^0®  =  3.2  MPa. 
The  dashed  line  simulation  uses  the  same  moduli  as  the  solid  line  simulation,  but  a  lower  maximum 
stress  of  (Tmax  =  2.8  MPa  and  less  severe  softening  parameters.  Note  that  all  simulations  exhibit 
a  peak  stress  of  3.2  MPa,  although  pre-peak  hardening  and  more  gradual  post-peak  softening  are 
observed  in  the  dashed  line  simulation. 

Although  pre-peak  hcirdening  can  be  modeled  with  the  present  implementation  of  the  brittle  dam¬ 
age  model,  as  just  demonstrated,  numerous  iterations  may  be  required  to  fit  the  model  to  data. 
This  is  because  the  user  specifies  a  threshold  of  cTmax  =  2.8  MPa  while  trying  to  simulate  a  peak 
stress  of  3.2  MPa.  In  addition,  the  model  lacks  flexibility  because  there  is  a  trade-off  between 


34 


Figure  3.2.  Effect  of  parameter  selection  on  hardening  and  softening  response. 

modeling  pre-peak  damage  and  the  severity  of  the  post-peak  softening  response.  To  overcome  these 
drawbacks,  the  modification  described  in  the  next  paragraph  is  proposed,  and  recommended  for 
future  implementation. 

In  our  proposed  modification,  the  user  specifies  two  additioned  par£imeters,  F  and  do.  The  parameter 
F  is  the  pre-peak  damage  threshold,  specified  as  a  fraction  of  the  maximum  stress  <Tmax,  where 
0  <  F  <  1.  The  parameter  do  is  the  amount  of  pre-peak  damage,  specified  as  a  frax:tion  of  the  total 
damage.  For  example,  F  =  0.6  and  do  =  0.15  indicate  that  damage  initiates  at  a  stress  of  0.6amax 
and  the  moduli  axe  degraded  15%  prior  to  attaining  a  peeik  stress  of  <Tmax- 


Modeling  pre-peak  damage  using  such  an  approach  requires  implementation  of  a  pre-peak  dam¬ 
age  function,  as  well  as  modification  of  the  damage  function  previously  given  in  Equation  2.5  for 
modeling  post-peak  response.  The  pre-peak  damage  function  is: 


d  =  ;?  Z?-do  for  Fro<n<  ro 


(3.4) 


(1  -  F)ro 

Note  that  d  in  Equation  3.4  varies  linearly  with  7^,  t.e.  d  =  0  at  =  Fro  and  d  =  do  at  fi  =  Tq. 
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Post-peaJc  degradation  is  readily  modeled  as  follows: 


i  =  do  +  (1  -  <i.)(l  -  <1^1^  -  Aexp®^*''”  - 

u 
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SECTION  4 

MODELING  CRACK  OPENING  AND  CLOSING 


We  also  implemented  a  feature  to  allow  recovery  (‘healing’)  of  the  bulk  modulus  during  cyclic 
loading  between  the  brittle  and  ductile  regimes.  This  feature  is  intended  to  simulate  opening  (bulk 
modulus  degradation)  and  closing  (bulk  modulus  recovery)  of  cracks.  Our  assumption  is  that  any 
existing  microcracks  will  close  under  compressive  pressure  (volumetric  strain),  but  will  initiate, 
grow,  and  coalesce  under  tensile  pressure  (volumetric  strain). 

Bulk  modulus  recovery  is  modeled  with  both  the  stress-based  and  strain-based  brittle  formulations 
previously  discussed  in  Sections  3.1  and  3.2.  For  the  stress-based  (strain-based)  model,  we  set 
the  brittle  damage  parameter  equal  to  zero  (‘healing’)  any  time  the  pressure  (volumetric  strain) 
becomes  compressive.  We  set  the  brittle  damage  parameter  back  to  it  previous  value,  md  allow 
further  brittle  damage  accumulation,  any  time  the  pressure  (volumetric  strain)  becomes  tensile. 

A  Roimd  D  benchmark  calculation  predicting  ‘healing’  during  dynamic  bending  of  a  reinforced 
concrete  slab  is  given  in  Section  8.  A  simpler  example  calculation  showing  cyclic  loading  between 
the  brittle  and  ductile  regimes  is  presented  in  Figure  4.1.  Loading  into  the  brittle  regime  and 
subsequent  degradation  are  simulated,  with  isotropic  damage  accumulation  of  d  =  0.3.  Moduli  K 
and  G  are  degraded  30%  from  their  initial  undamaged  values.  Now  continue  this  simulation  with 
loading  into  the  ductile  regime.  The  shear  modulus  G  will  retain  its  degraded  value,  but  the  bulk 
modulus  K  will  recovery  its  initial  und2unaged  value.  Hence  the  damage  is  no  longer  isotropic. 
Although  additional  degradation  of  G  may  occur  in  the  ductile  regime,  no  degradation  of  K  is 
allowed.  If  we  continue  this  simulation  with  loading  back  into  the  brittle  regime,  isotropic  damage 
accumulation  will  be  recovered.  For  this  example,  the  values  of  K  and  G  will  be  30%  less  than  their 
initial  undamaged  vedues. 
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SECTION  5 

MODELING  DILATION  DAMAGE 


Volume  expansion  under  compressive  loading  is  commonly  measured  in  concrete  specimens  during 
TXC  tests  with  little  or  no  confining  pressure.  Such  tests  also  exhibit  more  severe  softening  than  TXC 
tests  conducted  at  higher  confining  pressure.  Hence  it  is  possible  that  the  presence  of  volumetric 
tensile  strain  (or  possibly  lateral  tensile  strain)  contributes  to  the  more  severe  softening  observed 
in  TXC  tests  conducted  at  low  confining  pressure.  We  also  believe  that  dilation  behavior  is  critical 
in  Mode  II  (shear)  failure. 

We  implemented  a  dilation  damage  formulation  to  model  an  increase  in  the  severity  of  damage  when 
volume  expansion  occurs.  This  model  is  used  in  conjunction  with  the  stress-based  brittle  model 
previously  described  in  Section  3.1.  This  model  is  not  used  with  the  strain-based  brittle  model, 
because  the  strain-based  model  already  accounts  for  dilation  damage. 

Dilation  damage  initiates  when  the  pressure  is  compressive  (P  >  0),  the  volumetric  strain  is  ten¬ 
sile,  and  the  volumetric  strain  exceeds  a  dilation  damage  threshold.  The  damage  threshold  is 
^0  =  \jKtl  where  c„  max  =  <^maxlK  is  the  maximum  volumetric  strain  allowed  during  isotropic 
tension  simulations.  Two  softening  parameters  are  user  specified  {A  =  Ad  and  B  =  Bd  in  Equa¬ 
tion  2.5)  to  model  the  strain-softening  response.  The  effective  strain,  f  =  tj  in  Equation  2.5,  is 
fd  =  \jKtl  where  is  the  current  value  of  the  volumetric  strain.  Dilation  damage  can  be  neglected 
by  appropriate  selection  of  the  softening  parameters  {Ad  =  1  and  Bd  =  0). 

The  softening  responses  of  two  unconfined  compression  simulations  are  compared  in  Figure  5.1  with 
and  without  dilation  damage.  Note  that  more  degradation  is  modeled  with  dilation  damage  than 
without  dilation  damage. 
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Figure  5.1.  Compeirison  of  unconfined  compression  simulations  with  and  without  modeling  dilation 
damage. 
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SECTION  6 

ROUND  A  BENCHMARKS:  LABORATORY  TEST 

SIMULATIONS 


The  Round  A  calculations  consists  of  two  groups  of  calculations,  called  Groups  1  and  2.  The  Group 
1  calculations  are  single  element  fits  to  axial  stress-strain  data  from  TXC  and  TXE  tests  of  plain 
concrete,  for  confining  pressures  between  0  and  100  MPa.  These  calculations  demonstrate  our  fitting 
procedure  for  the  isotropic  damage  model.  The  Group  2  calculations  are  multi-element,  laboratory 
test  simulations  of  an  unconfined  TXC  test.  These  studies  use  the  fit  obtained  from  the  Group 
1  studies  and  demonstrate  that  rate-independent  multi-element  simulations  predict  more  severe 
softening  than  single  element  models.  Multi-element  simulations  also  predict  diagonal  damage 
patterns  and  splitting,  which  are  typical  failure  modes  of  concrete  test  specimens.  These  modes  are 
not  predicted  with  single  element  models.  Although  it  is  eflScient  to  fit  concrete  data  with  single 
element  models,  such  fits  should  be  considered  preliminary.  The  goodness  of  these  preliminary  fits 
should  be  evaluated  using  multi-element  lab  test  simulations. 

6.1  SINGLE  ELEMENT  FITS  TO  PIECEWISE-LINEAR  DATA 
REPRESENTATIONS. 

RDA  provided  the  benchmark  calculators  with  piecewise-linear  representations  of  concrete  laboratory 
test  data,  rather  than  with  the  actual  data  curves.  RDA’s  representations  of  the  data  will  be  referred 
to  as  ‘data’  through-out  the  remainder  of  this  section. 

The  single  element  in  the  Group  1  benchmark  calculations  is  modeled  10  cm  high  by  5  cm  square. 
The  element  was  loaded  first  in  hydrostatic  compression.  Then  the  element  was  further  compressed 
in  the  axial  direction  using  prescribed  axial  velocities,  while  the  confining  pressure  load  was  held 
constant.  The  calculations  were  run  at  a  slow  strain  rate  of  0.001/s  to  avoid  dynamic  effects. 

Single  element  fits  of  the  model  to  TXC  data  are  given  in  Figures  6. 1-6.6  and  to  TXE  data  in 
Figure  6.7.  The  fitting  procedure  is  described  in  detail  in  Appendix  D.  Four  main  features  of  the 
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Figure  6.1.  Fit  of  model  to  TXC  data  without  confining  pressure, 
fit  are  discussed  here: 

•  The  initial  elastic  response  following  hydrostatic  compression; 

•  The  ‘hardening’  behavior,  which  we  define  as  the  nonlinear  response  prior  to  the  peak  stress; 

•  The  peak  stress  and  the  strain  at  which  the  peaik  stress  occurs; 

•  The  ‘softening’  behavior,  which  is  a  decrease  in  stress  beyond  the  peak  stress  with  an  increase 
in  strain. 

The  data  in  these  figures  is  plotted  as  axial  stress  versus  axial  strain.  The  data  for  the  confined  TXC 
tests  is  trilinear  prior  to  the  peak  stress.  The  first  linear  segment  (not  apparent  at  low  confining 
pressure)  corresponds  to  the  intial  hydrostatic  compression  of  the  specimen.  The  second  linear 
segment  is  the  initial  elastic  response,  while  the  third  segment  provides  the  hardening  response  to 
peak  stress.  The  softening  behavior  following  the  peak  stress  is  bilinear  until  a  constant  residual 
stress  is  attained.  The  reader  should  keep  in  mind  that  ‘real  data’  does  not  behave  in  a  piece-wise 
linear  manner. 
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Figure  6.2.  Fit  of  model  to  TXC  data  for  5  MPa  confining  pressure. 
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strain  (Percent) 

Figure  6.4.  Fit  of  model  to  TXC  data  for  20  MPa  confining  pressure, 


Strain  (Percent) 

Figure  6.5.  Fit  of  model  to  TXC  data  for  50  MPa  confining  pressure. 
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RDA  provided  till  calculators  with  the  initial  elastic  moduli,  so  the  measured  and  calculated  responses 
are  in  good  agreement  during  the  initial  hydrostatic  compression  and  elastic  regimes,  although  dif¬ 
ferences  between  the  measured  and  calculated  responses  are  apparent  at  higher  confining  pressures. 
These  differences  are  due  to  the  calculation  of  plasticity,  which  initiates  at  70  MPa  for  hydrostatic 
compression,  or  slight  damage  accumulation. 

The  peak  stress  of  the  calculated  response  is  in  reasonable  agreement  with  the  data  for  most  fits. 
However,  the  calculated  strain  at  which  the  peak  stress  occurs  is  less  than  that  of  the  data,  for 
all  fits.  The  hardening  behavior  is  not  adequately  modeled.  This  is  probably  due  to  constraints 
imposed  by  the  dictation  of  elastic  moduli  and  initiation  of  plastic  response  as  implied  by  the 
piecewise-lineax  representations  of  data.  Good  fits  to  the  hardening  behavior  of  wSMR-5  concrete 
were  obtained  in  the  Round  B  benchmark  calculations. 

In  fitting  the  model  to  the  softening  behavior  of  the  data,  one  must  keep  in  mind  that  the  observed 
softening  in  TXC  laboratory  specimens  is  always  accompanied  by  localization  of  damage.  Test  spec¬ 
imens  typically  split  or  fail  on  a  diagonal,  and  the  distribution  of  strain  is  not  uniform  throughout 
the  specimen.  We  were  provided  with  the  average  axial  stress-strain  response  of  the  specimen,  which 
is  different  than  the  stress-strain  response  in  the  localized  damaged  regions  and  in  the  undamaged 
regions.  These  features  of  concrete  behavior  are  thoroughly  discussed  for  the  Group  2  laboratory 
test  simulations  in  Section  6.2. 

In  assigning  the  benchmark  calculations,  RDA  suggested  that  the  predicted  softening  in  the  Group  2 
unconfined  lab  test  simulations  will  be  steeper  than  the  Group  1  single  element  fits.  This  suggestion 
was  the  basis  of  our  Group  1  fit  for  unconfined  compression  test,  as  shown  in  Figure  6.1.  The 
softening  behavior  of  the  single  element  fit  is  not  nearly  as  steep  as  that  of  the  data.  However, 
for  confining  pressures  between  5  and  50  MPa  in  Figures  6.2-6. 5,  the  softening  behavior  of  the  fit 
is  similar  to  that  of  the  data.  For  a  confining  pressure  of  100  MPa,  the  fit  softens  while  the  data 
remains  perfectly  plastic  beyond  about  6.3%  strain. 

One  way  to  check  our  fit  of  the  damage  model  to  the  TXC  data  is  to  preform  three-dimensional  lab 
test  simulations  using  the  model  parameters  obtained  from  the  single  element  fit.  In  this  way  the 
average  axial  stress-strain  response  from  the  multi-element  lab  test  simulations  can  be  correlated 
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with  the  data.  One  such  correlation  is  made  for  the  unconfined  compression  test  and  discussed  in 
Section  6.3. 

We  also  fit  the  model  to  data  from  a  TXE  test  without  confinement,  as  shown  in  Figure  6.7.  The 
TXE  data  provided  by  RDA  is  a  tensile  strength  of  2.5  MPa  with  no  residual  strength  after  failure. 
At  the  time  these  calculations  were  performed,  a  tensile  cutoff  model  was  implemented  in  the 
concrete  damage  model.  With  the  tensile  cutoff  pressure  is  reached,  the  damage  is  instantaneous 
and  complete  t.c.  the  damage  parameter  is  set  equal  to  d  =  0.99. 

6.2  MULTI-ELEMENT  SIMULATIONS  OF  TXC  TESTS. 

The  Group  2  calculations  are  multi-element  simulations  of  an  unconfined  TXC  test.  They  demon¬ 
strate  the  effect  of  mesh  refinement  on  the  axial  stress-strain  response,  particularly  the  softening 
response,  and  on  the  localization  of  damage.  These  studies  use  the  fit  obtained  from  the  Group  1 
calculations  to  model  concrete  behavior.  First,  we  compare  the  average  axial  stress  versus  strain 
prediction  with  the  measured  data.  Next,  contours  of  damage  are  presented  and  indicate  how 
damage  localizes  into  diagonal  patterns.  Finally,  local  stress  versus  strain  predictions  are  examined 
inside  and  outside  the  damage  zones. 

The  specimens  are  modeled  10  cm  in  the  axial  z-direction,  and  5x5  square  cm  in  the  xy  plane. 
The  calculation  matrix  shown  in  Table  6.1  provides  the  mesh  refinement,  the  y-face  and  lateral 
edge  constraints,  and  whether  or  not  the  specimen  is  ‘pinched’  to  4.8  cm  in  the  x-direction  at 
midheight.  The  edge  constraints  of  all  pinched  specimens  are  free,  which  is  called  ‘lubricated’  in 
Table  6.1,  while  the  edge  constraints  of  all  unpinched  specimens  are  fixed.  For  both  the  pinched  and 
unpinched  specimens,  the  y-faces  are  either  free  or  in  plane  strain.  In  all  cases,  only  one  element  is 
modeled  in  the  y-direction.  We  will  refer  to  this  direction  as  the  ‘through-the-thickness’  direction. 

Three  mesh  refinements  are  also  used  in  the  Group  2  studies,  which  we  will  refer  to  as  coarse  for 
a  1x1x2  mesh,  medium  for  a  1x1x32  mesh,  and  fine  for  a  16x1x32  mesh.  The  sample  is  loaded 
by  a  compressive  axial  velocity  which  ramps  up  to  0.0015  m/ms  at  0.4  ms  and  remains  constant 
until  the  calculation  ends  at  1.5  ms.  The  final  axial  strain  is  1.95%.  Thus  the  average  strain  rate 
is  approximately  13/s,  which  is  much  higher  than  the  Group  1  strain  rate.  These  calculations 
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Table  6.1.  Matrix  of  Group  2  laboratory  test  simulations  for  Round  A. 


Case 

Zoning 

xxyxz 

Constraint 
on  y-Faces 

Lateral 
Constraint 
on  Ends 

Pinching 

free 

lubricated 

yes 

free 

fixed 

no 

1x1x2 

plane  strain 

lubricated 

yes 

2-CPF 

1x1x2 

plane  strain 

fixed 

no 

2-MUL 

free 

lubricated 

yes 

2-MUF 

free 

fixed 

no 

2-MPL 

1x1x32 

plane  strain 

lubricated 

yes 

2-MPF 

1x1x32 

plane  strain 

fixed 

no 

2-FUL 

16x1x32 

free 

lubricated 

yes 

2-FUF 

16x1x32 

free 

fixed 

no 

2-FPL 

16x1x32 

plane  strain 

lubricated 

yes 

2-FPF 

16x1x32 

plane  strain 

fixed 

no 

demonstrate  the  effect  on  mesh  refinement  on  the  average  axial  response  and  on  the  localization  of 
damage. 

Average  Axial  Stress-Strain  Response.  The  average  axial  stress-strain  response  from  the 
coarse,  medium,  and  fine  mesh  simulations  are  shown  in  Figures  6.8-6.10,  respectively.  Each  figure 
separately  shows  results  for  specimens  with  and  without  pinching.  For  the  fine  mesh  calculations, 
the  average  stress  is  calculated  by  averaging  the  axial  stress  for  all  16  elements  along  the  loaded 
edge  of  the  specimen.  The  average  strain  is  calculated  by  dividing  the  axial  displacement  along  the 
loaded  edge  by  the  10  cm  length  of  the  specimen. 

Plane  strain  and  unconfined  (free)  responses  are  given  in  each  plot  and  indicate  that,  in  all  cases, 
the  pealc  stress  is  higher  in  plane  strain  than  without  confinement.  This  is  because  the  peak  stress 
increases  with  confinement  in  TXC  tests,  so  the  model  provides  for  increasing  strength  with  confining 
pressure.  The  plane  strain  condition  provides  some  confinement,  i.e  the  through- the-thickness  stress 
is  not  zero,  so  the  peak  stress  is  enhanced. 

The  calculated  responses  of  the  pinched  and  unpinched  specimens  are  similar  for  the  medium  and 
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a)  Pinched  specimens  without  lateral  end  constraints. 


b)  Unpinched  specimens  with  lateral  end  constraints. 
Figure  6.8.  Coarse  mesh  computational  results. 
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a)  Pinched  specimens  without  lateral  end  constraints. 


b)  Unpinched  specimens  with  lateral  end  constraints. 
Figure  6.9.  Medium  mesh  computational  results. 
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Pinched  specimens  without  lateral  end  constraints. 


b)  Unpinched  specimens  with  lateral  end  constraints. 
Figure  6.10.  Fine  mesh  computational  results. 
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fine  mesh  studies,  but  not  for  the  coaxse  mesh  studies.  For  example,  the  fine  mesh  results  for  the 
pinched  specimens  in  Figure  6.10a  are  similar  to  those  of  the  unpinched  specimens  in  Figure  6.10b. 
However,  the  coeirse  mesh  results  for  the  pinched  specimens  in  Figure  6.8a  are  different  than  those 
of  the  impinched  specimens  in  Figure  6.8b.  Two  differences  are  noted:  the  magnitude  of  the  peah 
stress  and  the  oscillatory  behavior.  Both  of  these  differences  are  due  to  the  lateral  edge  constraint 
placed  on  the  impinched  specimens,  but  not  on  the  pinched  specimens. 

The  magnitude  of  the  peak  stress  in  the  coMse  mesh  studies  is  discussed  first.  All  nodes  except 
the  mid-height  nodes  are  constrained  from  latereJ  motion  in  the  unpinched  specimens  with  coarse 
mesh  refinement  (2  elements).  This  lateral  edge  constraint  provides  confinement,  which  increases 
the  peak  stress  in  the  unpinched  specimens.  The  peak  stress  is  greater  with  lateral  confinement 
(unpinched)  than  without  confinement  (pinched).  The  effect  of  the  lateral  edge  constraint  is  less 
significant  for  the  medium  and  fine  mesh  refinements,  because  only  two  of  the  33  nodes  along  each 
edge  of  the  mesh  are  laterally  constrained. 

Oscillations  are  also  apparent  in  the  pinched  specimen,  coarse  mesh  calculations,  but  not  in  the 
unpinched  specimen  calculations.  The  oscillations  in  the  coarse  mesh,  pinched  specimen  calculations 
can  be  eliminated  by  either  reducing  the  strain  rate,  or  by  applying  lateral  edge  constraints.  The 
axial  stress  oscillations  are  accompanied  by  synchronized  lateral  oscillations  in  both  stress  and 
displacement.  Without  lateral  end  constraints,  oscillatory  lateral  stresses  develop  because  the  lateral 
displacements  oscillate  about  their  quasi-static  values  at  the  Group-2  strain  rate.  Any  overshoot 
in  lateral  stress  is  accompanied  by  an  overshoot  in  axial  stress,  similar  to  the  effect  of  confining 
pressure  on  the  peak  axial  stress. 

One  point  to  note  about  the  the  computed  responses  is  that  the  steepness  of  the  softening  response 
increases  dramatically  with  mesh  refinement.  With  coarse  mesh  refinement,  the  average  axial  stress 
is  more  than  half  of  the  peak  stress  at  an  axial  strain  of  1%,  in  all  cases.  With  fine  mesh  refinement, 
the  average  axial  stress  is  less  than  10%  of  the  peak  stress  at  1%  strain,  in  all  cases. 

To  determine  how  well  we  fit  the  damage  model  parameters  to  the  data,  the  average  stress-strain 
prediction  from  the  the  unconfined  pinched  specimen  calculation  (2-FUL)  is  compared  with  the  test 
data  and  single  element  fit  in  Figure  6.11.  The  softening  response  of  the  lab  test  simulation  is 
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Figure  6.11.  Comparison  of  the  fine  mesh,  unconfined  lab  test  simulation  and  the  single  element  fit 
with  test  data. 

much  steeper  than  the  single  element  fit,  and  in  better  agreement  with  the  data.  This  comparison 
indicates  that  a  single  element  cannot  be  used  to  fit  the  softening  response  measured  from  a  lab  test 
simulation,  at  least  for  a  rate-independent  model,  Xhis  is  because  the  post-peak  strain  distribution 
is  highly  non-uniform  throughout  the  specimen,  due  to  localized  damage  accumulation,  and  this 
damage  accumulation  cannot  be  captured  with  a  single  element  model.  Damage  localization  is 
demonstrated  and  the  effect  of  further  mesh  refinement  is  discussed  in  subsequent  paragraphs. 

Another  difference  between  the  lab  test  simulation  and  single  element  fit  is  the  peak  stress  prediction. 
The  higher  peak  stress  calculated  with  the  lab  test  simulation  is  probably  due  to  dynamic  overshoot 
due  to  the  high  strain  rate.  The  lab  test  simulation  was  run  at  a  dynamic  strain  rate  of  13/s  whereas 
the  single  element  fit  was  run  at  a  much  slower  strain  rate  of  0.001 /s. 

Contours  of  Damage.  Contours  of  the  damage  parameter  d  are  plotted  in  Figures  6.12-6.15 
for  the  fine  mesh  calculations.  The  time  of  each  plot  is  after  substantial  post-peak  softening  has 
occurred.  A  value  of  d  =  0.99  indicates  that  meiximum  damage  has  occurred  and  the  moduli  have 
been  reduced  to  1%  of  their  original  elastic  values.  A  value  of  d  =  0  indicates  that  no  damage  has 
occurred. 
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contour  ualues 
R-  2.30E-01 
B-  3.8BE-01 
C-  4.89E-01 
D-  6.19E-01 
E-  7.4BE-01 
F-  B.7BE-01 


in-  r  1  cllsp.  scale  factor  -  0.1B0E+01  (ciefault)  , 

figure  0.12.  Contours  of  damage  for  an  unconnned  specimen  pinched  at  midheight  (2-FUL). 
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FPL  2 

time  -  0.59939E-«-e0 

contours  of  eff.  plastic  strain 


contour  ualues 


p- 

2.27E-01 
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3.46E-01 
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8.20E-01 
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2-FUr  MODEL  FOR  GROUP  2  BENChlM 
time  >  0.60000E+0Q 

contours  of  ef f .  plastic  strain 

min"  1.27SE—01  in  element  4 

max"  9.900E— 01  in  element  505 


contour  value 
0-  2.00E-01 
B=  2.90E-01 
C=  3.B0E-01 
D=  4.69E-01 
E"  5.59E-01 
F-  6.49E-01 
G>  7.3BE-01 
H"  0.2BE-01 
I"  g.lSE-01 


disp.  scale  factor 


0.100E+01  (default) 


Figure  6.14.  Contours  of  damage 
(2-FUF). 


for  an  unconfined  specimen  with  fixed  lateral  edge  constraints 
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2-FPF  MODEL  FOR  GROUP  2  BENCHM 
time  0.G00eBE-fe0 
contours  of  ef f .  plastic  strain 

min-  1.363E-01  in  element  6 

max-  9.900E-01  in  element  303 


contour  value 
A-  2.00E-01 
B-  2.97E-01 
C-  3.B6E-01 
D-  4.74E-01 
E-  5.63E-01 
F-  6.52E-01 
G-  7.41E-01 
H-  8.30E--01 
I-  g.lBE-01 


dlsp.  scale  factor  -  0.100E+01  (default) 


Figure  6.15.  Contours  of  damage  for  a  plane  strain  specimen  with  fixed  lateral  edge  constraints 

(2-FPF). 
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Figure  6.16.  The  deformed  configuration  of  a  plane  strain  specimen  pinched  at  mid-height  (2-FPL). 

The  damage  pattern  for  the  pinched  specimen  is  a  diamond  pattern.  Damage  initiates  at  the  mid¬ 
height  pinch  on  each  side,  and  at  center  on  the  loaded  and  unloaded  ends.  The  damage  pattern 
for  the  unpinched  specimens  with  lateral  end  constraints  is  an  ‘X’  pattern.  Damage  initiates  at  the 
four  laterally  constrained  corners  of  the  specimen. 

The  damage  pattern  can  also  be  detected  from  the  deformed  configuration  of  each  mesh,  although 
less  readily  than  from  the  contour  plots.  The  deformed  configuration  for  the  pinched  specimen 
without  confinement  is  shown  in  Figure  6.16.  Half  of  the  diamond  pattern  is  clearly  visible  on  the 
loaded  end  of  the  specimen.  Here  some  elements  have  distorted,  probably  because  their  moduli  have 
been  reduced  to  1%  of  their  original  values.  A  viscous  method  for  hourglass  stabilization  method 
was  used  in  these  analyses.  We  switched  to  a  stiffness  method  in  Round  B,  which  reduces  some  of 
the  observed  distortion. 

Local  Stress-Strain  Responses.  Here  we  examine  stress  and  strain  histories  at  specific  locations 
inside  and  outside  the  damage  zones.  Four  sets  of  stress  and  strain  histories  are  plotted  in  Figure  6.17 
for  the  plane  strain  specimen  pinched  at  mid-height  (2-FPL).  Three  of  these  locations  axe  at  the 
corner  points  of  the  diamond  pattern  region  of  severe  damage,  labeled  top  center,  side  pinch,  and 
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Figure  6.18.  Stress-strain  responses  at  four  different  locations  (2-FPL). 

bottom  center.  The  top  is  the  loaded  edge  of  the  mesh,  and  the  bottom  is  the  unloaded  edge.  The 
fourth  location  is  the  top  corner,  or  loaded  edge  comer. 

The  stress  histories  at  all  four  locations  are  similar,  indicating  that  the  distribution  of  axial  stress  is 
fairly  uniform  throughout  the  specimen,  even  while  damage  accumulates.  The  strain  histories  are 
similax  up  to  a  time  of  about  0.35  msec,  just  prior  to  the  time  of  peak  stress,  but  then  they  diverge. 
Axial  strains  as  large  as  13%  (center)  and  less  than  0.1%  (corner)  occur  on  the  loaded  edge  of  the 
specimen.  Thus  the  post-peak  distribution  of  strain  is  highly  non-uniform,  and  this  affects  the  local 
stress-strain  responses. 

Stress-strain  curves  are  plotted  in  Figure  6.18  at  the  same  four  locations.  The  stress-strain  curve  at 
top  center  reaches  a  peak  and  then  continuously  softens  as  the  strain  increases.  The  stress-strain 
curves  at  the  side  pinch  and  bottom  center  also  soften  following  the  peak  stress,  but  then  they 
unload  elastically  at  about  3%  strain.  The  stress-strain  curve  outside  the  severe  damage  region 
(top  corner)  loads  to  a  peak  stress  of  about  50  MPa  and  then  unloads  elastically.  This  comparison 
demonstrates  how  different  stress-strain  behaviors  are  predicted  at  various  locations  within  the 
specimen. 
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Time  (msec) 

a)  Stress  histories. 


Time  (msec) 

b)  Strain  histories. 

Figure  6.17.  Stress  and  strain  histories  calculated  at  four  different  locations  (2-FPL). 
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Figure  6.19.  Average  axial  response  calculated  with  two  different  tensile  cutoff  behaviors  (2-FUL). 

Effect  of  Tensile  Cutoff  Model  on  Calculated  Response.  The  average  axial  stress-strain 
curves  for  the  fine  mesh  calculations  in  Figure  6.10  exhibit  high  frequency  oscillations  after  substan¬ 
tial  softening  has  occurred.  This  high  frequency  behavior  occurs  when  the  tensile  cutoff  is  reached. 
We  implemented  the  cutoff  behavior  to  instantaneously  calculate  maximum  damage  (d  =  0.99) 
when  the  stress  state  reaches  the  tensile  cutoff.  However,  we  allowed  the  damage  parameter  to 
revert  to  its  previous  value  if  the  stress  state  subsequently  moves  off  of  the  cutoff  surface.  This 
resulted  in  abrupt  fluxuation  of  the  damage  parameter  which  caused  the  high  frequency  behavior. 

We  ran  an  additional  calculation  in  which  maximum  damage  (d  =  0.99)  was  sustained  for  the 
remainder  of  the  calculation  once  the  stress  state  reached  the  tensile  cutoff.  We  will  refer  to 
this  implementation  as  the  modified  tensile  cutoff  behavior.  Average  axial  stress-strain  histories 
calculated  with  the  original  and  modified  cutoff  behaviors  are  compared  in  Figure  6.19  for  the 
unconfined,  pinched  cylinder.  Note  that  most  of  the  high  frequency  behavior  is  eliminated  in  the 
calculation  with  the  modified  tensile  cutoff,  otherwise  the  stress-strain  curves  are  nearly  identical. 

Contours  of  damage  from  the  modified  tensile  cutoff  calculation  are  given  in  Figure  6.20  for  compar¬ 
ison  with  those  from  the  original  calculation  previously  shown  in  Figure  6.12.  Although  a  diamond 
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max-  0.990E+B0  in  element  511 
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P-  2.30E-01 
6-  3.60E-01 
C-  4.09E-01 
D-  6.19E-01 
E-  7.4BE-01 
F-  B.7aE-01 


Figure  6.20.  Contours  of  damage^  cSculated  with  the  modified  ten^le  cuto^  mo^el  ^2-FUL) 


pattern  of  damage  is  calculated  with  both  the  original  and  modified  tensile  cutoff  behaviors,  an 
axial  fracture  down  the  center  of  the  specimen  is  also  predicted  with  the  modified  behavior. 


6.3  EFFECT  OF  MESH  REFINEMENT  ON  THE  CALCULATED 
RESPONSE. 

Adding  Through-the-Thickness  Elements.  We  performed  a  three-dimensional  lab  test  simu¬ 
lation  to  determine  how  well  our  single  element  fit  predicts  the  average  axial  stress-strain  response 
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Figure  6.21.  Effect  of  through-the-thickness  refinement  on  the  average  axial  stress-strain  response 
of  an  unconfined  specimen  pinched  at  midheight. 

of  a  lab  test  specimen  in  unconfined  compression.  This  calculation  is  different  than  the  Group  2 
simulations  in  two  respects.  First,  the  calculation  was  run  at  a  slow  strciin  rate  of  0.01/s  to  eliminate 
dynamic  effects.  This  allows  us  to  more  readily  compare  the  lab  test  simulation  with  the  single 
element  fit  because  both  calculations  are  run  at  very  low  strain  rates.  Second,  elements  are  modeled 
through-the-thickness  of  the  specimen,  so  the  resulting  mesh  is  three-dimensional  and  should  be  a 
better  simulation  of  an  actual  test  than  the  Group  2  models  with  one  element  through-the-thickness. 
The  drawback  of  adding  elements  through-the-thickness  is  that  the  computer  run-time  is  increased. 


The  average  axial  stress-strain  response  is  shown  in  Figure  6.21  for  a  pinched  specimen  18  elements 
high  and  9x9  elements  over  the  square  cross  section.  This  calculation  is  referred  to  as  the  ‘3D 
simulation’.  The  average  stress  is  obtained  by  averaging  the  axial  stress  of  the  81  elements  across 
the  loaded  edge  of  the  specimen.  Note  that  the  peak  stress  and  pre-peak  response  calculated  with 
the  3D  mesh  are  nearly  equal  to  that  of  the  single  element  fit.  However,  the  post-peak  softening 
response  of  the  two  calculations  differs  dramatically.  The  softening  response  of  the  3D  calculation 
is  in  good  agreement  with  the  data,  whereas  the  softening  response  of  the  single  element  fit  is  much 
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less  steep  than  that  of  the  data.  These  comparison  suggest  that  a  rate-independent  fit  based  on  a 
single  element  model  will  not  adequately  capture  the  softening  response  of  a  lab  test  specimen. 

Also  shown  in  Figure  6.21  is  the  Group  2  simulation  with  one  element  modeled  through-the- 
thickness.  The  width  of  the  stress-strain  pulse  for  the  3D  simulation  is  less  than  that  for  the  Group 
2  simulation,  and  in  better  agreement  with  the  data.  This  may  be  because  through-the-thickness 
elements  were  modeled,  or  because  the  strain  rate  was  reduced  in  the  3D  simulation. 

Contours  of  damage  on  each  face  of  the  3D  mesh  are  shown  in  Figure  6.22  at  a  time  when  substantial 
softening  has  occurred  and  the  average  axial  strain  is  0.32%.  Note  that  although  the  contours 
indicate  a  diagonal  damage  pattern,  the  pattern  is  not  a  diamond  pattern,  nor  is  it  symmetric 
about  the  axial  midplane  of  each  face.  The  reason  for  the  lack  of  symmetry  is  not  known,  although 
it  may  be  do  to  sensitivity  of  the  damage  and  tensile  cutoff  portions  of  the  model  to  numerical 
roimdoff. 

Refining  the  Group  2  Fine  Mesh.  To  assess  the  effect  of  mesh  refinement  on  the  fine  mesh 
predictions,  an  additional  calculation  for  simulation  FUL-2  was  performed  by  doubling  the  number 
of  elements  in  each  coordinate  direction.  The  resulting  mesh  is  64  elements  high  and  32x2  elements 
over  the  cross  section.  Stress-strain  results  for  the  original  16x1x32  mesh,  and  the  refined  32x2x64 
mesh  are  compared  in  Figure  6.23  and  demonstrate  that  the  pre-pealc  and  post-peak  responses  are 
nearly  identical  until  the  concrete  softens  to  about  60%  of  its  peaJc  stress. 

Contours  of  damage  from  the  refined  mesh  calculation  are  shown  in  Figure  6.24  for  comparison  with 
those  from  simulation  FUL-2  and  previously  shown  in  Figure  6.12.  Again,  a  diamond  pattern  of 
damage  is  predicted,  although  two  possible  fracture  lines  are  calculated  with  the  refined  mesh  while 
one  possible  fracture  line  is  calculated  with  the  original  mesh.  However,  the  width  of  the  damage 
bands  calculated  with  the  refined  mesh  is  approximately  equal  to  that  of  the  original  mesh. 
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contour  voluos 
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Figure  6.22.  Contours  of  damage  on  each  face  of  the  3D  simulation. 


Figure  6.23.  Eifect  of  mesh  refinement  on  the  average  axial  stress-strain  response  of  an  unconfined 
specimen  pinched  at  midheight  (2-FUL). 
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SECTION  7 


ROUND  B  BENCHMARKS:  FITS  TO  WSMR-5  DATA 


The  Round  B  calculations  are  similar  to  the  Roimd  A  calculations  except  that  RDA  provided  the 
calculators  with  WSMR-5  data  curves  measured  by  WES.  The  Round  B  calculations  consists  of 
Group  1  single  element  fits  and  Group  2  multi-element  simulations.  However,  the  Round  B  calcu¬ 
lations  are  run  at  a  slower  strain  rate  than  the  Roimd  A  calculations  to  eliminate  the  undesirable 
oscillatory  behavior  observed  in  the  Round  A  simulations.  The  Round  B  calculations  also  model 
rate  dependence  of  the  strength  using  the  Duvaut-Lions  viscoplastic  update  algorithm  previously 
discussed  in  Section  2.4.4. 

7.1  SINGLE  ELEMENT  FITS  TO  WSMR-5  DATA. 

Single  element  fits  of  the  model  to  WSMR-5  data  are  given  in  Figures  7.1  through  7.5.  The  calcula¬ 
tions  were  run  at  a  strain  rate  of  0.001/s  to  eliminate  inertial  effects,  and  rate  dependence  of  the 
strength  was  neglected.  We  used  the  isotropic  damage  model  previously  described  in  Section  2.1. 
We  fit  the  data  without  considering  damage  first,  then  iteratively  adjusted  the  fits  to  include  dam¬ 
age,  as  discussed  in  Appendix  D.  Table  7.1  provides  the  material  model  parameters  used  in  these 
fits. 

We  obtained  excellent  fits  of  the  model  to  the  isotropic  compression  and  uniaxial  strain  data.  We 
also  modeled  less  softening  than  measured  from  the  TXC  tests.  This  is  because  we  expected  more 
softening  to  be  predicted  in  the  Group  2  TXC  simulations  than  modeled  with  the  single  element 
fits.  The  direct  pull  test  data  was  fit  to  the  tensile  pressure  cutoff  model  previously  described  for 
the  Round  A  benchmarks. 

A  three-dimensional  calculation  was  performed  to  asses  the  effect  of  mesh  refinement  on  the  cal¬ 
culated  response.  This  calculation  was  not  part  of  the  required  benchmark  series.  The  three- 
dimensional  simulation  was  performed  with  a  16x16x32  mesh  for  comparison  with  the  single  element 
fit,  as  shown  in  Figure  7.6.  Rate  dependence  of  the  strength  was  ignored.  The  softening  response 
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Table  7,1.  Model  parameters  for  the  fit  with  damage. 


Constitutive  Model 


K  (MPa) 

G  (MPa) 

p  (gm/cm^) 

2.30 

Failure  Envelope 


F’*(Ji)  =  a-i 
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Figure  7.1.  Fit  of  the  model  to  isotropic  compression  data. 


Figure  7.2.  Fit  of  the  model  to  uniaxial  strain  data. 
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Figure  7.3.  Fit  of  the  model  to  principal  stress  versus  pressure  data  for  a  uniaxial  strain  test. 

of  the  three-dimensional  simulation  degrades  more  rapidly  than  the  single  element  fit,  although  the 
degradation  is  not  ais  rapid  as  we  would  have  liked  to  have  modeled^.  Although  no  data  is  available 
on  the  softening  response  of  WSMR-5  concrete  from  the  unconfined  compression  tests,  WES’s  recom¬ 
mended  curve  softens  to  zero  strength  at  about  0.5%  strain.  These  comparisons  suggest  that  more 
damage  should  be  modeled  in  our  fit  to  the  WSMR-5  data.  They  also  suggest  that  the  softening 
response  is  more  sensitive  to  the  mesh  refinement  when  more  damage  is  modeled. 

7.2  MULTI-ELEMENT  SIMULATIONS  OF  TXC  TESTS. 

The  calculation  matrix  in  Table  7.2  provides  the  mesh  refinement  and  y-face  constraints  for  the 
Round  B  Group  2  calculations.  Only  the  lubricated  end  (free)  cases  from  Round  A  were  re¬ 
calculated  in  Round  B.  The  Round  B  multi-element  simulations  were  run  at  an  average  strain  rate 
of  c  =  1.3,  which  is  a  factor  of  ten  slower  than  the  Round  A  multi-element  simulations.  The  slower 
strain  rate  eliminates  undesirable  oscillations.  The  Round  B  simulations  were  also  performed  with 
strength  enhancement  as  a  function  of  strain  rate.  The  strength  enhancement  factor  previously 
given  in  E<juation  2.44  was  provided  by  RDA  to  all  benchmark  participants. 

^Time  constraints  prevented  us  from  re-adjusting  our  fits  and  performing  additional  simulations. 
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PRINCIPAL  STRESS  DIFFERENCE 


Figure  7.6.  Comparison  of  three-dimensional  simulation  with  the  single  element  fit  (Round  B). 


Table  7.2.  Matrix  of  Group  2  laboratory  test  simulations  for  Round  B. 


Case 

Zoning 

xxyxz 

Constraint 
on  y-Faces 

Lateral 

Constraint 

on  Ends 

Pinching 

B2-CUL 

1x1x2 

free 

lubricated 

yes 

B2-CPL 

1x1x2 

plane  strain 

lubricated 

yes 

B2-MUL 

1x1x32 

free 

lubricated 

yes 

B2-MPL 

1x1x32 

plane  strain 

lubricated 

yes 

B2-FUL 

16x1x32 

free 

lubricated 

yes 

B2-FPL 

16x1x32 

plane  strain 

lubricated 

yes 
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The  average  axial  stress-strain  responses  from  the  coarse,  medium,  and  fine  mesh  simulations  are 
compared  in  Figure  7.7a  for  plane  strain  and  in  Figure  7.7a  for  unconfined  compression.  The  peak 
stresses  attained  with  the  coarse,  medium,  and  fine  mesh  simulations  are  in  close,  but  not  identical 
agreement.  As  expected  from  Round  A,  the  peak  stress  attained  in  the  plane  strain  simulations  is 
greater  than  the  peak  stress  attained  from  the  unconfined  compression  simulations.  The  value  of 
the  peak  stress  attained  in  all  simulations  is  greater  than  the  static  value  of  45  MPa  because  rate 
dependence  of  the  strength  was  modeled. 

Unlike  the  Round  A  simulations,  the  steepness  of  the  softening  response  for  the  Round  B  simulations 
does  not  increase  dramatically  with  mesh  refinement.  Only  small  increases  in  softening  response 
are  calculated  with  increased  mesh  refinement,  at  least  for  plane  strain  response.  This  is  because 
we  modeled  less  damage  in  the  Round  B  simulations  than  in  the  Round  A  simulations.  This 
demonstrates  that  the  effect  of  mesh  refinement  on  the  softening  response  is  less  severe  when  less 
damage  is  modeled. 

RDA  compiled  all  participants  Round  B  benchmark  calculations  in  Reference  9.  APTEK  s  original 
B2-FUL  simulation  reported  to  RDA  exhibits  a  sudden  drop  in  stress.  This  instability  (sudden  stress 
drop)  is  associated  with  the  type  of  hourglass  control  we  used.  The  sudden  stress  drop  is  due  to 
instability  of  the  rectangular  mesh  caused  by  hourglassing^  in  elements  with  aspect  ratios  of  16. 
Such  a  high  aspect  ratio  is  outside  the  usual  range  of  2  or  less. 

Two  types  of  hourglass  stability  methods  are  available  in  DYNA3D,  called  viscous  and  stiffness 
methods.  The  default  form  in  DYNA3D  is  a  viscous  form.  The  effect  of  the  hourglass  stability 
method  on  simulation  B2-FUL  is  shown  in  Figure  7.8.  Note  that  a  sudden  drop  in  stress  occurs 
when  using  the  viscous  form  for  hourglass  stability,  but  not  when  using  the  stiffness  form.  Our 
original  (unstable)  B2-FUL  simulation  reported  to  RDA  in  Reference  9  was  preformed  with  the 
viscous  method  whereas  the  stable  response  was  preformed  with  the  stiffness  method.  All  other 
Round  A  and  B  calculations  were  performed  with  the  viscous  method. 


^Hourglass  deformation  modes  are  spurious  zero  energy  modes  that  arise  from  elements  formulated  with  one-point 
Gauss  quadrature.  The  spurious  deformation  modes  must  be  stabilized  while  retaining  legitimate  deformation  modes. 
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Axial  Strain  (Percent) 


(a)  Plane  strain. 


Axial  Strain  (Percent) 

(b)  Unconfined  compression. 

Figure  7.7.  Comparison  of  coarse,  medium,  and  fine  mesh  simulations  for  Round  B  lubricated  spec 


Figure  7.8.  Effect  of  hourglass  stability  method  on  the  unconfined,  fine  mesh  calculation  (B2-FUL). 
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SECTION  8 


ROUND  D  BENCHMARKS:  DYNAMICS  OF 
REINFORCED  SLABS 

The  Round  D  benchmarks  study  the  dynamics  of  reinforced  concrete  slabs  in  plane  strain  under 
various  load  levels.  The  geometry  of  the  undeformed  slab  mesh  in  shown  in  Figure  8.1.  The  bottom 
surface  of  the  slab  coincides  with  the  x  —  y  plane  with  plane  strain  in  the  y-direction.  The  vertical 
z-direction  corresponds  to  the  through-the-thickness  direction  of  a  concrete  wall.  The  mesh  is 
10x1x6  solid  concrete  elements  in  the  x-y-z  directions,  respectively.  Thus  6  elements  are  modeled 
through-the-slab  thickness.  The  dimensions  of  the  slab  are  380x25x114  mm.  The  bottom  and  top 
rows  of  elements  are  13  mm  high.  All  other  elements  are  22  mm  high. 


Traction  Free 

Figure  8.1.  Geometry  of  slab  in  Round  D  benchmark  calculations. 

Two  rebar  are  are  bonded  to  the  concrete,  one  between  the  first  and  second  rows  of  elements,  and 
one  between  the  fifth  and  sixth  rows  of  elements.  The  cross-sectional  area  of  each  rebar  is  17.2 
mm^.  The  rebar  are  located  on  the  front  face  at  j/  =  0. 

The  left  edge  at  r  =  0  mm  is  clamped.  The  right  edge  at  x  =  380  mm  is  a  roller  boundary.  It 


78 


Figure  8.2.  Unsealed  pressure  history  applied  to  slab  sections  in  Round  D. 
Table  8.1.  Case  designation  and  load  factor  for  Round  D  pressure  pulse. 


Case 

Load  Factor 
LF 

D1 

0.6 

D2 

1.4 

D3 

2.2 

D4 

3.0 

represents  the  midspan  symmetry  plane  of  a  concrete  wall.  The  front  and  back  surfaces  at  y  =  0 
mm  and  j/  =  25  mm  are  roller  boundaries  to  enforce  the  plane  strain  condition.  The  bottom  surface 
at  2  =  0  rnm  is  traction  free.  The  top  surface  at  2  =  114  mm  has  a  spatially  uniform  pressure 
history,  p{t),  applied.  The  pressure  history  consists  of  a  sharp  spike  followed  by  a  more  gradual 
pulse,  as  shown  in  Figure  8.2.  The  magnitude  of  this  pressure  history  is  scaled  by  the  load  factors 
given  in  Table  8.1  for  four  different  cases.  Before  being  scaled,  the  sharp  spike  initiates  at  0  ms, 
peaks  to  8.25  MPa  at  0.01  ms,  then  decays  to  0  Mpa  at  0.02  ms.  The  more  gradual  pulse  initiates 
at  0.79  ms,  peaks  to  0.6  Mpa  at  0.80  ms,  then  decays  to  0  MPa  in  21.6  ms. 

RDA  compiled  all  benchmark  participants  results  and  reported  them  in  Reference  11.  The  thrust, 
bending  moment,  shear  force,  and  deflection  histories  compiled  by  RDA  for  APTEK  axe  shown  in 
Figure  8.3.  The  four  calculations  were  run  for  20  msec.  Of  particular  interest  are  the  displacement 
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histories  calculated  along  the  midspan  symmetry  plane  at  x=380  mm  and  z  =  hi  mm.  The 
maximum  downward  displacement  ranges  from  0.026  cm  for  LF=0.6  to  13.9  cm  for  LF=3.0,  which 
is  almost  three  orders  of  magnitude  in  mciximum  displacement  for  a  factor  of  5  difference  in  scale 
factor.  The  time  to  reach  the  maximum  displacement  also  increases  with  load  factor. 

The  mayiTmiTn  displacement  versus  load  factor  and  time  to  maximum  displacement  are  compared 
in  Figure  8.4  for  all  participants  models.  These  comparisons  demonstrate  the  extreme  sensitivity 
of  the  maximum  deflection  to  small  changes  in  the  applied  load.  Most  calculations  predict  a  large 
jump  in  deflection  for  load  factors  between  2.2  and  3.0,  just  as  APTEK  does. 

Contours  of  damage  are  shown  in  Figure  8.5  through  Figure  8.8  for  cases  D1  through  D4,  respec¬ 
tively.  In  all  cases,  tensile  damage  initiates  on  the  clamped  edge  near  the  top  (loaded)  surface  and 
along  the  bottom  surface  near  the  midspan.  Very  limited  damage  is  predicted  in  calculation  D1  by 
20  ms.  The  location  of  this  damage  corresponds  to  the  front  face  (loaded  face)  of  a  concrete  wall 
near  the  edge.  Much  more  severe  damage  is  predicted  in  calculation  D2  by  20  ms,  particularly  over 
the  bottom  surface  of  the  slab.  This  surface  corresponds  to  the  backface  of  a  concrete  wall. 

Contours  of  damage  are  predicted  at  6  ms  and  20  ms  for  calculations  D3.  The  damage  contours 
predicted  at  20  ms  are  not  significantly  different  than  those  predicted  at  6  ms.  These  contours 
suggest  that  backface  spall  is  likely  to  occur  near  the  midspan. 

Contours  of  damage  at  6,  8.5,  11,  and  20  msecs  are  shown  for  calculation  D4.  The  entire  slab 
is  severely  damaged  by  8.5  ms.  Note,  however,  that  damage  accumulation  is  reduced  near  the 
clamped  edge  by  11  ms.  This  is  due  to  modeling  modulus  recovery  or  ‘healing’  of  the  bulk  modulus, 
as  previously  discussed  in  Section  4.  Initially,  elements  along  the  clamped  edge  go  into  tension 
(P  <  0)  and  both  the  bulk  and  shear  moduli  are  severely  degraded.  As  the  slab  deforms  further, 
some  of  the  elements  along  the  clamped  edge  go  into  compression  {P  >  0).  In  this  case,  the  bulk 
modulus  is  recovered,  or  set  back  to  its  initial  undamaged  value.  Although  this  example  illustrates 
our  modulus  recover  model,  the  model  still  must  be  verified  against  experimental  data. 
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Max  Displacement 


Figure  8.4.  Maximum  deflection  versus  load  factor,  compiled  by  RDA  for  Round  D  slab  sections 


Figure  8.5.  Contours  of  damage  at  20  ms  for  Round  D  case  D1  calculation. 


Figure  8.6.  Contours  of  damage  at  20  ms  for  Round  D  case  D2  calculation. 
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SECTION  9 


SUMMARY  AND  FUTURE  EFFORTS 


The  three-invariant  smooth-cap  model  with  isotropic  damage  models  the  essential  features  of  con¬ 
crete  behavior.  These  features  include  different  strengths  in  TXC  and  TXE,  shear  enhanced  com¬ 
paction,  dilatancy,  post-peak  softening,  degradation  of  the  elastic  moduli,  and  localized  damage 
accumulation.  The  damage  modeling  efforts  discussed  in  this  report  are  a  first  attempt  at  modeling 
softening  and  modulus  degradation,  and  all  results  should  be  considered  preliminary.  Numerous 
enhancements  and  studies  are  recommended  for  a  follow-on  effort,  as  follows: 

•  Stress-Based  vs.  Strain-Based  Models.  Simo  and  Ju  developed  alternative,  but  not 
equivalent,  damage  formulations  for  stress-based  and  strain-based  models.  Although  we  im¬ 
plemented  their  strain-based  formulation  for  ductile  response,  future  efforts  will  examine  the 
appropriateness  of  stress  vs.  strain-based  formulations,  particularly  for  modeling  cyclic  re¬ 
sponse  between  the  brittle  and  ductile  regimes.  Whichever  type  of  formulation  is  selected,  it 
will  be  used  consistently  for  both  brittle  and  ductile  damage. 

•  Continuity  Between  Brittle  and  Ductile  Models.  Currently  two  separate  models  are 
implemented  for  the  brittle  (tensile  pressure)  and  ductile  (compressive  pressure)  regimes. 
Although  the  damage  initiation  threshold  is  modeled  continuously  between  these  two  models, 
the  evolution  of  the  damage  parameter  and  the  accumulation  of  plasticity  are  discontinuous 
between  these  two  models  at  P  =  0.  Future  efforts  will  ensure  that  there  is  a  continuous  and 
smooth  connection  between  modeling  damage  in  the  brittle  and  ductile  regimes  by  merging  the 
brittle  and  ductile  models  into  a  single,  more  sophisticated  model.  This  model  will  contain 
a  pressure-dependent  damage  initiation  threshold  throughout  the  compressive  and  tensile 
regimes,  allow  for  slight  plasticity  in  the  tensile  regime  (as  warranted  by  available  data),  and 
smoothly  adjust  the  severity  of  the  softening  response  (damage  accumulation)  between  the 
tensile  and  compressive  regimes. 
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•  Isotropic  vs.  Anisotropic  Damage  Models.  In  addition  to  the  isotropic  damage  formula¬ 
tion  discussed  in  this  report,  Simo  and  Ju  also  developed  an  anisotropic  damage  formulation 
for  modeling  the  tensile  failure  of  concrete.  Unlike  the  scaler  damage  variable  used  in  the 
isotropic  model,  a  fourth-order  tensor  characterizes  damage  in  the  anisotropic  model.  Future 
efforts  will  examine  the  feasibility  of  implementing  their  anisotropic  model  into  our  concrete 
damage  model.  We  also  plan  to  review  the  brittle  anisotropic  model  recently  implemented 
into  DYNASD  by  Govindjee  in  Reference  20. 

•  Modeling  the  Effects  of  Preloading.  The  measured  reduction  in  modulus  and  strength 
due  to  preloading  a  specimen  in  uniaxial  strain  prior  to  direct  pull  or  unconfined  compression 
tests  was  was  previously  shown  in  Figure  l.*4  for  WSMR-5  concrete.  Future  efforts  will  predict 
the  measured  stress-strain  responses  and  damage  mechanisms.  The  main  difficulty  here  is 
determining  the  evolution  of  a  pressure-dependent  damage  threshold  (through  the  softening 
parameters  A  and  B)  when  cycling  between  tension  and  compression.  In  addition,  the  use 
of  stress-based  or  strain-based  damage  models  will  have  a  significant  affect  on  the  calculated 
results. 
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APPENDIX  A 

THREE  INVARIANT  SCALING  FUNCTION  TZ 


Rubin  considers  a  three-invariant  yield  function  of  the  general  form: 


(A.l) 


where  /  =  0  at  fjulure  and  Yd  is  the 
Rubin  is 


failure  surface  for  TXC.  The  functional  form  of  Tl  proposed  by 


11  = 


-61  +  ^6? -46260 
262 


(A.2) 


where 


62  =  (cos  ^  —  a  sin +  6sin^ 

61  =  a(cos  —  a  sin  /3) 

60  =  -(3 -I- 6  -  o")/4 
b  =  (2Qi  +  o)^  —  3 

,  (A.3) 

— Ol  +  VOl  —  40200 

02  =  Q2 

oi  =  +  2Qi(Q2  —  1) 

ao  =  2Ql{Q2-l) 


The  value  of  H  in  Equation  A.2  depends  on  the  state  of  stress  through  the  angle  and  on  exper¬ 
imentally  determined  values  for  Qi  and  Q2  as  functions  of  pressure.  The  function  H  provides  the 
shape  of  the  failure  envelope  in  the  octahedral  plane.  This  shape  may  be  that  of  a  circle  (moSS), 
semi-regular  hexagon  (Mohr- Coulomb),  or  irregular  hexagon-like  shape  in  which  each  of  the  six 
sides  is  quadratic  (rather  than  linear)  between  the  TXC  and  TXE  stress  states. 
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The  yield  function  previously  given  in  Equation  2.18  for  the  APTEK  concrete  model  is  of  the  general 
form  suggested  by  Rubin  in  Equation  A.l,  where  Yd  =  V^Fj  for  Fc  =  1.  The  general  form  of  K 
given  in  Equation  A. 2  has  been  implemented  into  DYNaSd  for  pressure  dependent  values  oi  Qi  and 
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A.l  REDUCTION  OF  11  TO  THE  MOHR-COULOMB  CRITERION. 

Mohr’s  condition  for  pressure-dependent  materials  states  that  yielding  occurs  when  the  maximum 
shear  stress  reaches  a  critical  value  as  a  function  of  pressure,  see  Reference  15.  Rubin  derives  an 
analytical  expression  for  It  for  a  nonlinear  failure  envelope  as  follows: 


W)  = 


_ ^3^2 _ 

(1  +  Q2)  cos  ^  -  \/3(l  -  Qi)  sin 


(A.4) 


where 


Qi 

Q2 


V^Q2 

1  +  <52 

3  —  sin  ip 
3  -h  sin 


(A.5) 

(A.6) 

(A.7) 


and  the  friction  angle  ip  may  be  pressure  dependent. 

The  simplest  form  of  Mohr’s  condition  is  known  as  the  Mohr-Coulomb  criterion  and  is  graphically 
represented  in  Figure  A.l  as  a  straight-line  envelope  of  fcdlure  states  in  normal  stress,  <t,  versus 
shear  stress,  t,  space.  The  principal  stresses  are  >  <72  >  <73.  Two  experimentally  determined 
constants  define  the  straight-line  envelope:  the  cohesion,  c,  and  the  friction  angle,  (^,  as  follows: 


|t|  =  c  —  crtany? 


(A.8) 


The  principal  stresses,  and  hence  lr|,  can  be  defined  in  terms  of  the  three  stress  invariants,  Ji,  yjj‘2, 
and  The  failure  envelope  also  forms  a  straight  line  in  vs  Ji  stress  invariant  space. 
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Figure  A.l.  Mohr-Coulomb  criterion  with  a  straight-line  failure  envelope. 

APTEK’s  three-invariant  yield  function,  Equation  2.18,  is  equal  to  a  straight-line  Mohr-Coulomb 
criterion  if  TZ  is  given  by  Equation  A.4,  and  the  failure  surface  parameters  a  and  6  are  related  to 
the  friction  angle,  (f,  and  unconfined  compressive  strength,  (T„,  as  follows: 


6c  cos 

■\/3(3  —  sin^) 

0  _  2  sin  V? 

\/3(3  —  sinv?) 

where 

"  2VN 

^  ^  1  +  siny> 

1  —  sin 


(A.9) 

(A.IO) 


(A.11) 

(A.12) 


The  failure  surface  parameters  and  TZ  depend  on  the  same  pressure-independent  friction  angle.  For 
most  soils,  the  friction  angle  is  y  «  30  degrees. 


One  point  to  note  is  that  the  three-invariant  model  is  more  general  than  a  strict  Mohr-Coulomb 
criterion.  As  previously  noted  for  the  Mohr-Coulomb  criterion,  the  shape  of  the  failure  envelope 
in  the  J\  plane,  and  the  shape  of  the  semi-regular  hexagon  in  the  octahedral  plane,  depend 
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on  the  same  friction  angle  <p.  ip  is  constant  for  a  streiight-line  envelope  and  a  function  of  pressure 
for  a  non-linear  envelope.  In  the  three-invaricint  model,  p  may  be  used  to  determine  the  shape  of 
the  irregular  hexagon  in  the  octahedral  plane.  The  failure  surface  function,  Fj,  may  be  specified 
independently  of  p,  and  does  not  have  to  be  linear  even  though  a  constant  value  for  p  is  specified. 

A.2  REDUCTION  OF  Tl  TO  THE 

MAXIMUM-OCTAHEDRAL-SHEAR  STRESS  CRITERION. 

This  criterion  states  that  yielding  begins  when  the  octahedral-shear  stress  reaches  a  critical  value. 
Rubin  assumes  a  slight  generalization  of  the  MOSS  condition  by  allowing  the  critical  value  to  be  a 
function  of  pressure.  Rubin  demonstrates  that,  for  this  case,  the  function  Tl  reduces  to  7^  =  (5i  = 
^2  =  1"  Thus  the  failure  curve  Yd  is  not  scaled  for  different  states  of  stress.  The  yield  function  in 
Equation  A.l  reduces  to  a  two-invariant  yield  function.  For  this  generalized  criterion,  the  failure 
envelope  is  nonlinear  with  pressure,  but  remains  a  circle  in  the  octahedral  plane. 

The  Drucker-Prager  and  Von  Mises  criteria  are  special  cases  of  the  generalized  MOSS  criterion.  The 
straight-line  envelope  of  the  Drucker-Prager  criterion  forms  a  cone  in  principal  stress  space,  while 
the  pressure-independent  Von  Mises  envelope  forms  a  cylinder. 

A.3  REDUCTION  OF  TO  A  WILLAM-WARNKE  MODEL. 

The  Wiliam- Warnke  model  is  a  subset  of  the  more  flexible  Rubin  model  previously  discussed  in 
Section  2.4.2.  In  these  models,  the  shape  of  the  failure  surface  in  the  deviatoric  plane  transitions 
with  pressure.  This  means  that  the  ratio  of  the  TXE  to  TXC  failure  surface  increases  with  pressure, 
as  does  the  ratio  of  the  TOR  to  TXE  failure  surfaces.  In  a  Wiliam- Warnke  model,  the  TOR/txe 
ratio  depends  specifically  on  the  TXE/txc  ratio  at  each  pressure  value.  In  the  more  general  Rubin 
model,  the  TOR/txe  ratio  can  be  specified  independently  of  the  TXE/txc  ratio. 

The  Wiliam- Warnke  model  was  the  specification  for  the  Round  D  and  E  benchmarks.  However,  no 
TXE  or  TOR  data  was  measured  for  WSMR-5  concrete.  Instead  of  fitting  our  model  to  data,  RDA 
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specified  the  location  and  shape  of  the  TXE  surface  along  the  hydrostat,  denoted  $2,  as  follows: 


0.686  P  <  75  MPa 

^2{P)  =  <  0.610  +  0.00102P  75  MPa  <  P  <  382  MPa  (A.13) 

1.0  P  >  382  MPa 

In  addition,  RDA  specified  that  the  shape  of  the  failure  surface  in  the  deviatoric  plane  must  transition 
with  pressure  according  to  a  Willam-Wamke  model. 

To  model  the  Willam-Wamke  failure  surface,  we  fit  the  functional  form  of  our  TXE  surface,  Q2(P), 
to  RDa’s  specification  in  Equation  A.13.  Then  we  derived  the  location  of  the  Willam-Wamke  TOR 
surface,  denoted  'i'l,  from  the  following  formula: 


^i(P)  = 


\/3(l  —  Ql)  +  {2Q2  —  1)  *  ^3(1  —  -F  5(^2  —  4(^2 
3(1 -<31)  +  (1- 2^2)2 


(A.14) 


where  Q2  is  a  function  of  pressure.  This  formula  can  readily  be  reduced  from  textbook  formulas  for 
radius  of  the  Wiliam- Warnke  model  in  the  deviatoric  plane  (see,  for  example.  Reference  15).  Finally 
we  fit  the  functional  form  of  our  TOR  surface,  Qi{P),  to  the  derived  specification  in  Equation  A.14. 


The  Rubin  model  can  accommodate  a  wide  range  of  three-invariant  model  behaviors.  Although 
our  model  meets  the  Wiliam- Warnke  benchmark  specifications,  we  have  the  flexibility  to  adjust  our 
model  if  those  specifications  turn  out  to  be  inappropriate  for  high-strength  concrete.  For  example, 
if  TXE  and  TOR  concrete  data  is  measured,  our  model  can  be  directly  fit  to  that  data,  even  if  the 
data  does  not  transition  according  to  a  Willam-Wamke  model. 
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APPENDIX  B 

EXPLICIT  PORE  COMPACTION  MODEL  FOR  SOILS 


In  the  CONWEB  test  series  in  Reference  21,  high  stress  levels  are  generated  in  soil  in  the  vicinity 
of  an  explosive  charge.  Mechanical  property  data  for  backfill  materials  used  in  these  tests,  i.e., 
dynamic  uniaxial  stress-strain  data  indicate  that  the  backfill  materials  stiffen  at  high  stress  levels. 
Hence  they  are  not  adequately  modeled  with  a  constant  bulk  modulus  model.  In  consultations 
with  Professor  Miles  Rubin  of  the  Israel  Technion  Institute,  it  was  decided  that  a  Mie-Gruneisen 
equation  of  state,  coupled  with  explicit  treatment  of  pore  collapse  at  low  pressure,  would  be  an 
appropriate  representation  for  both  low  and  high  pressure  states.  The  Mie-Gruneisen  equation 
and  fitting  procedures  to  Hugoniot  shock  data  are  well  established.  We  have  integrated  it  into 
the  smooth-cap  model  to  treat  in  a  continuous  manner  both  high  and  low  pressure  behavior.  The 
theory  of  the  Mie-Gruneisen  model  and  evaluation  of  the  implementation  are  discussed  in  the  next 
two  subsections. 

B.l  THEORY. 

Consider  a  geologic  material,  such  as  soil,  which  consists  of  solid  particles  with  voids  between 
particles,  called  ‘pores’.  Let  V  be  the  initial  volume  of  the  soil  in  the  unstressed  state,  and  v  the 
current  volume  in  the  stressed  state.  The  total  volume  of  the  soil  is  the  sum  of  the  volumes  of  the 
solid  particles,  denoted  by  the  subscript  ‘s’,  and  the  pores,  denoted  by  the  subscript  ‘p’,  as  follows; 

V  = 

V  =  Va  +  Vp  (B.l) 

The  relative  volume  of  the  soil,  J,  is  the  ratio  of  the  current  to  initial  volume,  J  =  v/V. 

The  pressure  in  the  soil  is  taken  to  depend  on  the  pressure  in  the  solid  particles,  Pa,  and  the  soil 


B-1 


porosity,  <^,  in  the  following  form: 


P  =  [1-<^(«)]P.(£;„J,«)  (B.2) 

in  which  k  is  the  cap  hardening  parameter  and  Eg  is  the  specific  internal  energy  per  unit  mass 
in  the  solid  particles,  calculated  from  the  strain  energy  and  mass  of  the  material.  Heat  transfer 
is  neglected.  The  soil  porosity  is  the  ratio  of  the  current  pore  volume  to  the  current  soil  volume, 
<f>  =  Up/u.  The  initial  porosity  of  the  soil  is  $  =  Vp/V.  Note  that  the  pressure  in  the  soil  is  equal 
to  the  pressure  in  the  solid  particles  when  the  porosity  is  zero,  i.e.,  P  =  Ps  when  =  0.  Until  the 
pores  are  fully  collapsed,  the  soil  pressure  is  less  than  the  solid  particle  pressure,  i.e.,  P  <  F,  when 
0. 

To  determine  the  pressure  in  the  soil,  expressions  for  the  solid  particle  pressure  and  the  porosity 
must  be  defined.  The  pressure  in  the  solid  particles  is  given  by  the  Mie-Gruneisen  equation  of 
state.  The  porosity  is  defined  in  terms  of  the  cap  hardening  relation.  This  relation  describes  the 
plastic  volume  strain  as  a  function  of  «.  Expressions  for  these  quantities  are  given  the  following 
two  sections. 

B.1.1  Mie-Gruneisen  Equation  of  State. 

The  Mie-Gruneisen  equation  gives  the  pressure  in  a  material  as  a  function  of  internal  energy  and 
density,  related  to  pressures  and  energies  along  a  reference  curve  at  the  same  density.  The  form  of 
the  Mie-Gruneisen  equation  given  in  Reference  22  is: 


P-Pref  =  /)r(P-P,ef)  (B.3) 

where  P  is  pressure,  E  is  internal  energy,  p  is  density,  and  the  subscript  ‘ref’  indicates  the  reference 
state,  r  is  a  material  parameter  called  the  Gruneisen  ratio,  and  is  obtained  from  laboratory  test 
data.  It  is  is  assumed  to  be  a  function  of  the  density  only;  in  particular,  pT  =  constant  in  the 
present  theory. 

The  Mie-Griineisen  equation  extends  information  from  known  pressure-volume  and  energy-volume 
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states,  typically  the  Hugoniot  states  for  shocked  material,  to  other  values  of  pressure  and  energy. 
With  the  Hugoniot  for  the  solid  particles  as  the  reference  state,  denoted  by  the  subscript  ‘hs’,  the 
pressure  in  the  solid  particles  is  given  by: 


The  pressure  and  energy  in  the  solid  particles  along  the  Hugoniot  are  given  by: 

POb^^cs 


Phs  = 


(1  -  s,<t>cBy 


2(1  -  si<i>cBy 


(B.5) 

(B.6) 


where  po,  is  the  initial  density  of  the  soUd  particles,  co  is  the  bulk  sound  speed  in  the  solid  at  low 
pressures,  and  Si  is  the  slope  of  a  linear  fit  to  shock  velocity  versus  particle  velocity  data  from 
shock  tests.  The  constant  product  p^V,  is  equal  to  the  product  of  the  initial  values,  posBo,.  The 
parameter  <I>cb  is  the  relative  change  in  the  volume  of  the  solid  particles,  given  by: 


<l>ca 


1- 


Izi 

i_$ 


j 


(B.7) 


Note  from  Equations  B.5  and  B.6  that  the  Hugoniot  pressure  and  energy  in  the  solid  particles  are 
zero  when  the  change  in  volume  of  the  solid  particles  is  zero. 


B.1.2  Soil  Porosity. 

An  expression  is  derived  in  this  section  for  <^(k),  which  is  the  soil  porosity  as  a  function  of  the  cap 
hardening  parameter.  This  is  accomplished  in  two  steps.  First,  an  expression  is  given  for  the  soil 
porosity  as  a  function  of  the  the  plastic  volume  strain.  Second,  a  cap  hardening  relation  is  defined 
which  provides  the  plastic  volume  strain  as  a  function  of  k. 

It  is  assumed  that  the  plastic  volume  strain  of  the  soil  is  a  function  of  the  porosity  only  in  the 
following  form: 
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The  initial  porosity  of  the  soil  in  the  unstressed  state  is  $,  and  the  initial  plastic  volume  strain 
is  zero.  As  the  soil  is  compressed,  the  pores  collapse  until  the  porosity  is  zero.  At  this  state,  the 
plastic  volume  strain  is  maximum  and  is  equcd  to  the  initial  porosity.  No  further  plastic  volume 
strain  is  allowed.  Thus,  from  Equation  B.8, 

fP  =  0  when  4>  =  ^ 

=  $  when  4>  =  0  (B-9) 

Rearrangement  of  Equation  B.8  gives  the  current  porosity  as  a  function  of  the  plastic  volume  strain; 

(B.IO) 

1  —  Cv 

The  cap  hardening  equation  provides  a  relation  between  the  plastic  volume  strain,  ej,  and  the 
hardening  parameter,  k.  For  the  smooth-cap  model  this  relation  is: 

el  =  VF{i  -  (B  11) 

The  parameter  PF  =  $  is  the  maximum  plastic  volume  strain.  The  function  is  the  value  of 

the  first  stress  invariant,  Ji,  at  the  intersection  of  the  cap  with  the  Ji  axis.  The  parameter  Xq  is 
the  value  of  X  evaluated  at  /cq,  the  initial  value  of  the  hardening  parameter.  The  parameters  D 
and  D2  are  obtained  by  fitting  Equation  B.ll  to  experimental  data. 

B.2  EVALUATION  OF  THE  IMPLEMENTATION. 

Stress-strain  histories  calculated  with  the  Mie-Gruneisen  model  are  shown  in  Figure  B.l.  These 
comparisons  are  for  McCormick  Ranch  Sand  during  loading  and  unloading  in  uniaxial  strain.  Prop¬ 
erties  of  McCormick  Ranch  Sand  reported  in  Reference  3,  and  values  assumed  for  Fos  and  Sj,  are 
given  in  Table  B.l. 

The  response  of  the  Mie-Gruneisen  model  during  loading  and  unloading  in  the  low  pressure  region 
is  compared  in  Figure  B.la  with  the  response  calculated  with  a  constant  bulk  modulus  model.  The 
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strain 

(a)  Low  pressure  behavior. 


Strain 


(b)  High  pressure  behavior. 

Figure  B.l.  Stress-strain  histories  calculated  with  Mie-Gruneisen  model  implemented  in  DYNA3D 
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Table  B.l.  Parameter  values  for  Mie-Gruneisen  demonstration  problem. 


Fo{Ji 


Constant  Bulk  Modulus  Model 


K  (ksi) 

66.7 

G  (ksi) 

40.0 

p  (lbf/in3) 

0.07 

Mie-Gruneisen  Model 


Co  (in/s) 

19250 

To, 

2.0 

Si 

1.5 

G  (ksi) 

40.0 

p  (Ibf/in^) 

0.07 

Failure  Envelope 


F,W  =  a  - 

-f  0  Ji 

a  (psi) 

250 

e 

0 

A  (psi) 

180 

^  (psi“^) 

0.67  X  10*3 

Moveable  Cap 
(Ji-I)(|Jx-Z|-f  Ji-X) 
■  2{X  -  Lf 


X(«)  =  X  -f  RFe{L) 
L{k)  =  max(Ko, «) 


R 

2.5 

Xo  (ksi) 

0.0 

Cap  Hardening  Rule 
;P  =  VF{1  - 


W 

0.066 

D  (psi“^) 

0.67  X  10-3 

£»2  (psi“^) 

0.0 
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responses  axe  in  good  agreement,  which  validates  the  low  pressure  behavior  of  the  Mie-Griineisen- 
smooth-cap  implementation. 

During  compressive  loading  to  about  10%  strain,  coUapse  of  the  pores  and  compression  of  solid 
particles  occurs  when  the  active  portion  of  the  model  is  the  cap  mode,  designated  by  ‘C  in  Fig¬ 
ure  B.la.  Collapse  of  the  pores  reduces  the  porosity  of  the  soil.  All  plastic  volume  strain  is  due  to 
the  change  in  porosity.  In  this  example  problem,  all  pores  are  completely  collapsed,  the  porosity  is 
zero,  and  the  plastic  volume  strain  is  maximum  (ej  =  W)  at  about  10%  strain.  The  partitioning  of 
total  strain  into  elastic  and  plastic  components  depends  on  user  selection  of  the  parameters  of  the 
cap  and  failure  surfaces. 

During  loading  or  unloading  in  the  elastic  regime,  designated  by  ‘E’  during  unloading  in  Figure  B.la, 
changes  in  pore  and  solid  particle  volume  also  occur,  but  in  such  a  manner  that  the  porosity  remains 
constant.  Hence  no  changes  in  plastic  volume  strain  occur. 

The  response  of  the  Mie-Gruneisen  model  in  the  high-pressure  region  is  shown  in  Figure  B.lb. 
Loading  and  unloading  occur  along  the  same  stress-strain  curve  for  strains  greater  than  about  10%, 
i.e.,  all  incremental  strains  are  elastic  and  due  to  deformation  of  the  solid  particles.  Below  10% 
strain,  loading  and  unloading  are  along  different  stress-strain  paths,  similar  to  the  low  pressure 
behavior  previously  shown  in  Figure  B.la. 
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Figure  C.l.  Axisymmetric  compression  of  a  Mohr-Coulomb  medium  around  a  circular  hole. 

APPENDIX  C 

VERIFICATION  OF  THE  THREE-INVARIANT  MODEL 

Kennedy  and  Lindberg  23  and  Florence  and  Schwer  24  provide  the  stress-strain  response  of  a  Mohr- 
Coulomb  medium  for  axisymmetric  compression  around  a  circular  hole.  The  medium  is  subjected 
to  an  increasing  free-held  pressure,  Fo>  and  a  constant  internal  pressure,  /*,■,  as  shown  schematically 
in  Figure  C.l.  A  plane  strain  condition  is  imposed.  Yielding  initiates  on  the  inner  boundary  of  the 
hole.  The  radius  Rp  of  the  yielded  region  increases  with  increasing  free-held  pressure. 

Analytical  results  are  calculated  for  a  straight-line  (constant  friction  angle)  yield  criterion.  Material 
properties  used  in  the  analysis  are  given  in  Table  C.l.  Analytical  stress-load  and  strain-load  histories 
axe  shown  in  Figure  C.2.  The  response  is  elastic-perfectly  plastic,  with  the  plastic  response  initiating 
at  about  Po  =  11  ksi.  These  histories  are  calculated  at  a  radius  of  r  =  2  inches  from  the  center  of 
the  hole.  The  value  of  the  normalized  invariant,  J3,  and  hence  11,  varies  throughout  the  elastic  and 
plastic  regimes. 

A  constant  bulk  modulus,  Kq,  is  used  to  calculated  the  pressure-volumetric  strain  response  of 
the  medium.  In  addition,  the  elastic  volumetric  strain,  is  the  difference  between  the  total 
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(a)  Stress-load  histories. 


(a)  Strain-load  histories. 

Figure  C.2.  Analytical  results  for  the  axisymmetric  compression  of  a  Mohr-Coulomb  medium  around 
a  circular  hole. 
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Table  C.l.  Material  properties  used  in  analysis  of  a  Mohr-Coulomb  medium. 


V 

0.3 

G  (ksi) 

1000 

(psi) 

2000 

<f>  (degrees) 

30 

Pi  (psi) 

500 

a  (inch) 

1.0 

r  (inch) 

2.0 

volumetric  strain,  ekk,  and  the  plastic  volumetric  strain,  elk  (the  double  subscript  indicates  the  use 
of  the  summation  convention).  The  analytical  pressure- volumetric  strain  response  is  as  follows. 


P  =  Koiekk  -  4k)  (^*1) 

This  is  a  different  pressure  relation  than  used  in  APTEK’s  explicit  pore  collapse  model,  so  the  analyt¬ 
ical  results  cannot  be  used  to  directly  verify  the  implementation  of  the  third  invariant  in  the  explicit 
pore  collapse  model.  However,  the  analytical  results  can  be  used  to  verify  the  implementation  of 
the  third  invariant  in  the  constant  bulk  modulus  model. 

Comparisons  Between  the  Analytical  and  Constant  Bulk  Modulus  Models.  In  the 
DYNa3d  calculations,  the  analytical  strain  histories  shown  in  Figure  C.2(b)  are  applied  to  a  single 
element^.  DYNA3d  output  stress  histories  axe  compared  with  the  analytical  stress  histories  shown 
previously  in  Figure  C.2(a).  Material  parameters  used  in  the  calculations  are  given  in  Table  C.2  and 
correspond  to  those  in  Table  C.l.  The  DYNa3d  failure  surface  parameters  a  and  0  are  related  to  the 
friction  angle,  v?,  and  unconfined  compressive  strength,  o-„,  used  in  the  analyses  (see  Equations  A.9 
and  A. 10). 

As  shown  in  Figure  C.3,  the  numerical  results  are  in  excellent  agreement  with  the  analytical  re¬ 
sults  throughout  the  elastic  and  plastic  regimes.  This  validates  the  extension  of  the  constant  bulk 
modulus  model  to  three  invariants.  In  particular,  it  validates  all  modifications  to  yield  function, 
return  mapping  algorithm,  and  the  stress  update.  These  same  modifications  are  implemented  in 

^They  are  interpreted  by  the  dyna3d  material  model  driver  as  strain-time  histories. 
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Table  C.2.  Material  properties  used  in  three-invariant  extensions  of  APTEKs  models. 


Bulk  Modrdus  Model 

Mie-Griineisen  Equation  of  State 

Ko  (ksi) 

2167 

Co  (in/s) 

109700 

(psi) 

0 

r 

0 

K2  (psi) 

0 

Si 

0 

G  (ksi) 

1000 

G  (ksi) 

1000 

P  (Ibf/in^) 

0.070 

p  (Ibf/in^) 

0.07 

Failure  Envelope 

a  (psi)  692.8 

0  0.231 

A  (psi)  0 

P  (psi"^)  0 

Cap  and  Hardening  Rule 

R  0.1 

Xq  (psi)  1  X  10^ 

$  0.05 

D  (psi-i)  0.67  X  10-^ 

£>2  (psi-^)  0.0 


0  5  10  15  20 


Load  (psi)  E3 

Figure  C.3.  Comparison  of  analytical  results  using  a  Mohr-Coulomb  criterion  with  numerical  results 
using  a  three-invariant  model  with  constant  bulk  modulus. 
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Figure  C.4.  Comparison  of  stress-load  histories  calculated  with  the  explicit  pore  collapse  and  con¬ 
stant  bulk  modulus  models  extended  to  three  invariants. 

the  explicit  pore  collapse  model. 

Comparisons  Between  the  Explicit  Pore  Collapse  and  Constant  Bulk  Modulus  Models. 
Comparisons  between  the  three-invariant  extensions  these  cap  models  are  given  in  Figure  C.4.  All 
properties  were  previously  given  in  Table  C.2.  Stress-load  histories  are  in  agreement  throughout 
the  elastic  response  and  during  the  early  portion  of  the  plastic  response.  However,  the  histories 
diverge  with  increasing  load  because  the  models  use  different  pressure  relations. 

The  pressure  relation  used  in  the  explicit  pore  collapse  model  is: 

P  =  (1  -  (0.2) 

where  <f>  is  the  porosity  and  P,  is  the  pressure  (refer  to  Appendix  B).  The  solid  particle  pressure 
is  represented  by  the  Mie-Gruneisen  equation  of  state.  Through  appropriate  selection  of  model 
parameters  (F  =  0,  5/  =  0),  Equation  C.2  is  reduced  to  a  model  with  a  constant  solid-particle  bulk 
modulus,  Kos- 

P  =  (1  -  <f>)Kos(l  -  Js)  (C.3) 

where  J$  is  the  relative  volume  of  the  solid  particles.  The  term  1  —  Js  is  the  change  in  volume  of 
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the  solid  particles  relative  to  the  initial  solid  particle  volume. 

There  are  three  differences  between  the  pore  collapse  and  constant  bulk  modulus  (analytical)  pres¬ 
sure  relations  given  by  Equations  C.3  and  C.l,  respectively: 

1.  The  models  have  different  bulk  moduli.  In  the  pore  collapse  model,  the  modulus  of  the  solid 
particles,  Kos,  is  constant.  The  effective  modulus,  K,  depends  on  the  porosity  as  well  as  the 
solid  particle  modulus,  namely:  K  =  {1  -  <f>)K».  Thus  the  effective  modulus  varies  between 
an  initial  value  of  (1  —  ^)K,  before  any  pores  are  compressed,  and  a  final  value  of  Kt  once 
all  pores  are  fully  crushed.  $  is  the  initial  porosity.  In  the  analytical  model,  the  effective 
modulus,  Kq,  is  constant.  It  is  related  to  the  solid-particle  modulus  and  initial  porosity  of 
the  pore  collapse  model,  as  follows:  Kq  =  (1  — 

2.  The  models  calculate  dilatation  differently.  In  the  pore  collapse  model,  the  total  relative 
volume,  t7,  is  accurately  calculated  from  the  determinant  of  the  deformation  gradient  matrix. 
In  turn,  the  solid  particle  volume,  J,,  depends  on  J  and  (j>  as  follows:  J,  =  J(1  —  ^)/(l  —  $).  If 
the  relative  volume  is  defined  as  J  =  1  —  A,  then  A  is  the  change  in  volume  per  unit  volume, 
and  is  called  the  dilatation  (see  Reference  25).  In  the  analytical  model,  the  total  volumetric 
strain  is  a  first  order,  small  strain  approximation  of  the  dilatation:  ejtjt  «  A. 

3.  The  models  have  different  definitions  of  plastic  volume  strain.  In  the  pore  collapse  model,  the 
plastic  volume  strain  is  the  difference  between  the  relative  volumes  of  the  solid  particles  and 
medium:  el  =  (J,  -  J)/Js-  In  the  analytical  model,  the  plastic  volume  strain  is  the  difference 
between  the  total  and  elastic  volume  strains:  ej  =  Ckk  ~ 

To  compare  the  pore  collapse  model  in  Equation  C.3  with  the  analytical  model  in  Equation  C.l,  each 
equation  is  rewritten  in  terms  of  ej,  K,  and  J.  The  constant  modulus  relation  from  Equation  C.l 
becomes: 

P  =  Ko{l-J-el)  (C.4) 

where  the  approximation  ejtfc  «  1  —  J  was  used.  The  pore  collapse  pressure  relation  from  Equa- 
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tion  C.3  becomes: 


(C.5) 


P  =  Ko 


l-J-tl 

(1  -  dy 


where 


!-<!> 


(C.6) 


These  relations  are  identical  when  the  plastic  volume  strain  is  zero,  such  as  during  the  initial  elastic 
response  (see  Figure  C.4),  but  diverge  as  the  plastic  volume  strain  increases  with  increasing  load. 
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APPENDIX  D 

FITTING  PROCEDURE 


This  section  describes  our  procedure  for  fitting  the  concrete  damage  model  to  data.  We  fit  the 
isotropic  damage  model  previously  discussed  in  Section  2.1  to  piecewise-Uneax  representations  of 
data  provided  by  RDA  for  the  Roimd  A  benchmarks.  The  actual  data  curves  were  not  provided. 
RDA’s  representations  of  the  data  will  be  referred  to  as  ‘data’  through-out  the  remainder  of  this 

section. 

First  we  describe  our  fit  of  the  model  to  data  without  including  damage.  Such  a  fit  predicts  elastic- 
perfectly  plastic  response  for  TXC  tests,  which  is  not  representative  of  the  data.  Next  we  describe 
our  fit  of  the  model  to  data  with  damage  included  in  the  fit. 

D.l  FIT  WITHOUT  DAMAGE. 

We  began  our  fitting  procedure  by  fitting  the  the  smooth-cap  model  to  data  without  considering 
damage.  The  parameters  we  obtained  from  this  fit  axe  given  in  Table  D.l.  Note  that  a  pressure 
cutoff  if  listed.  This  is  because  these  fits  were  made  prior  to  implementation  of  the  brittle  damage 
model.  Also  note  that  only  two  parameters  axe  listed  for  the  three  invariant  surface.  This  is  because 
a  reduced  form  of  the  Rubin  three  invariant  model  was  implemented  in  the  concrete  damage  model 
at  the  time  this  fits  were  made.  The  reduced  form  did  not  allow  the  shape  of  the  failure  surface  in 
the  deviatoric  plane  to  transition  with  pressure.  The  fitting  procedure  without  damage  is  standard 
so  only  a  brief  outline  is  given  here,  as  follows: 

•  The  initial  bulk  and  shear  moduli,  K  and  G  respectively,  were  provided  by  RDA.  In  general, 
however,  they  can  be  attained  from  unloading  curves  for  from  isotropic  compression  and 
uniaxial  strain  data  at  very  low  pressure; 

•  The  tensile  cutoff  was  obtained  from  the  unconfined  TXE  data.  This  fit  is  given  in  Figure  6.7; 
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Figure  D.2.  Preliminary  fits  without  damage  to  pressure-volumetric  strain  data. 

•  The  four  parameters  of  the  TXC  feiilure  surface  were  obtained  from  an  analytical  expression 
provided  by  RDA  for  the  peak  stress  difference  versus  pressure.  This  preliminary  fit  is  given 
in  Figure  D.la; 

•  The  two  parameters  for  the  three-invariant  model  were  obtained  from  an  analytical  expression 
for  the  TXE  surface,  and  a  straight-line  fit  between  the  TXE  and  TXC  values,  i.e.  straight-line 
Mohr-Coulomb  behavior  in  the  deviatoric  plane.  The  fit  to  TXE  data  is  given  in  Figure  D.lb; 

•  The  three  parameters  of  the  hardening  rule  were  obtained  from  the  elastic  bulk  modulus  and 
a  bilinear  pressure-volumetric  strain  curve  obtained  from  hydrostatic  compression  data.  This 
preliminary  fit  is  shown  in  Figure  D.2; 

•  The  initial  location  of  the  cap  was  obtained  from  the  isotropic  compression  pressure- volumetric 
strain  curve  while  the  cap  shape  function  was  obtained  from  a  fit  to  the  stress  path  from 
uniaxial  strain  data.  This  preliminary  fit  is  given  in  Figure  D.3. 

Results  of  the  fit  without  damage  are  compared  with  piecewise-linear  representations  of  TXC  test 
data  in  Figure  D.4.  The  confining  pressure  is  indicated  next  to  each  curve.  Although  we  fit  the  peak 
stress  fairly  well,  the  calculated  response  is  elastic-perfectly  plastic,  which  is  not  representative  of 
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Table  D,l.  Model  parameters  for  the  fit  without  damage. 


Constitutive  Model 


K  (MPa) 

14700 

G  (MPa) 

11025 

p  (gm/cm^) 

2.28 

Failure  Envelope 

F’e(Ji)  =  a-' 

a  (MPa) 

230 

e 

0 

A  (MPa) 

221 

/3  (MPa-i) 

1.6  X  10-3 

Moveable  Cap 

F(J  .A  1  + 

Fc{Ji,  K)-l  2(X  -  !)■•' 

A'(k)  =  L  +  RFe(L) 
L{k)  =  mox(«o,  k) 


R 

1.2 

Xo  (MPa) 

210 

Cap  Hardening  Rule 


W 

0.12 

D  (MPa-^) 

1.35  X  lO-'* 

Di  (MPa-2) 

1.00  X  10-® 

Third  Invariant 


Qx  =tor/txc  Q2  =txe/txc 


02 

0.705 

0.686 

Pressure  Cutoff 

T  (MPa) 

-2.5 

Damage  Model 


T 


Tq 

mg^ 

A 

■HH 

B 

HH 
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Figure  D.l.  Preliminary  fit  without  damage  to  TXC  and  TXE  failure  surfaces. 
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Figure  D.3.  Preliminary  fit  without  damage  to  the  stress  path  from  uniaxitil  strain  data. 

the  post-peak  softening  response  of  the  lab  test  data.  In  addition,  we  calculate  less  strain  hardening 
than  mecisured. 

D.2  FIT  WITH  DAMAGE. 

Our  fits  of  the  damage  model  to  TXC  data  are  given  in  Figure  D.5  for  confining  pressure  between  0 
and  100  MPa.  These  fits  should  be  compared  with  the  TXC  data  previously  given  in  Figure  D.4b. 
Note  that  the  model  predicts  the  softening  behavior  of  the  data  reasonably  well,  but  not  the 
hardening  response. 

The  fit  of  the  model,  with  and  without  damage,  to  the  unconfined  compression  data  is  shown  in 
Figure  D.6.  The  ratio  of  the  stress  from  the  fit  with  damage,  to  the  stress  from  the  fit  without 
damage,  is  the  reduction  factor  1  -  d.  The  point  at  which  the  fitted  curves  separate,  which  is  at 
about  25  MPa  in  this  example,  depends  on  the  value  of  Vg.  We  attempted  to  model  a  small  amount 
of  hardening  in  this  example  by  selecting  Vg  to  initiate  damage  prior  to  the  peak  stress.  We  could 
have  selected  to  initiate  damage  at  the  peak  stress.  However,  damage  would  still  have  initiated 
prior  to  the  peah  stress  in  the  fits  with  confining  pressure.  Thus  both  pre-peak  hardening  and 
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Stress  (MPa) 


Strain  (Percent) 


(a)  Fit  without  damage. 


Strain  (Percent) 

(b)  Data  provided  by  RDA. 

Figure  D.4.  Single  element  stress-strain  curves,  fit  without  damage  to  TXC  data. 


Figure  D.5.  Single  element  stress-strain  curves,  fit  with  damage  to  TXC  data. 

post-peak  softening  aie  modeled  with  the  damage  model,  although  the  hardening  response  is  not 
adequately  predicted  with  the  selected  damage  parameters. 

The  peak  stress  is  attained  by  iterative  adjustment  of  the  shear  failure  surface  location  and  damage 
model  parameters.  Two  fits  of  the  shear  failure  surface  to  the  principal  stress  difference  versus 
pressure  data  are  given  in  Figure  D.7.  The  fit  labeled  ‘fit  without  damage’  provides  the  peak  stress 
of  the  elastic-perfectly  plastic  curves  previously  given  in  Figure  D.4.  The  fit  labeled  ‘fit  with  damage’ 
provides  the  peak  stresses  calculated  with  and  without  damage  and  previously  shown  in  Figure  D.6. 
Note  that  we  adjusted  the  shear  failure  surface  upward  to  include  damage  in  the  model.  Otherwise, 
the  peak  stresses  attained  from  our  fit  with  damage  would  be  less  than  the  desired  values.  Our 
adjusted  parameters  for  the  fit  with  damage  axe  given  in  Table  D.2.  It  is  evident  from  our  fitting 
procedure  that  the  peak  stress  of  our  damage  model  depends  on  damage  accumulation  as  well  as 
plasticity. 

The  value  of  the  damage  parameter  d  depends  on  the  elastic  strain  energy  term,  t.  Equivalent 
strain  histories  are  shown  in  Figure  D.8  and  demonstrate  that  the  equivalent  strain  in  unconfined 
compression  is  four  times  larger  than  that  for  a  50  MPa  confining  pressure.  The  equivalent  strain 
decreases  with  confinement  in  TXC  calculations,  apparently  due  to  less  lateral  expansion  with  higher 
confinement,  so  the  damage  accumulation  also  decreases  with  confinement.  This  is  why  less  softening 
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Figure  D.6.  Fits  with  and  without  damage,  to  unconfined  compression  data. 


Figure  D.7.  Adjustment  of  the  shear  failure  surface  fit  to  account  for  pre-peak  damage. 
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Table  D.2.  Model  parameters  for  the  fit  with  damage. 


Constitutive  Model 


K  (MPa) 

14700 

G  (MPa) 

11025 

p  (gm/cm®) 

2.28 

Failure  Envelope 

Fe(Ji)  =  a- 

+  0Ji 

a  (MPa) 

12 

$ 

0.35 

A  (MPa) 

0 

^  (MPa-i) 

0 

FciJi, 


Moveable  Cap 

(Ji-L){\Ji-L\  +  J,-L) 

2{x  -  ly 


^(/c)  =  i  +  RFe{L) 
L(k)  =  max{Ko,  /c) 


R 

1.2 

Xo  (MPa) 

210 

Cap  Hardening  Rule 
eP  =  Vr{l  - 


W 

0.12 

D  (MPa-i) 

1.35  X  lO-'* 

Di  (MPa-2) 

1.00  X  10"® 

Third  Invariant 
Qi  =tor/txc  Q2  =txe/txc 


ai 

0.705 

OC2 

0.686 

Pressure  Cutoff 

T  (MPa) 

-2.5 

=  1  - 

Damage  Model 
ro(l  -  A)  ^^Biro  -  r 
r 

ro 

0.15 

A 

0.8 

B 

0.1 
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Figure  D.8.  Variation  of  equivalent  strain  energy  with  confining  pressure. 

is  observed  in  the  TXC  stress-strain  curves  as  the  confinement  level  increases. 

Two  fits  of  the  damage  model  to  the  unconfined  compression  data  are  shown  in  Figure  D.9.  The 
parameter  A  was  selected  to  fit  the  peak  stress  and  the  parameter  B  provides  the  steepness  of 
the  softening  response.  We  selected  Tq  to  initiate  damage  just  prior  to  the  peak  stress.  This  figure 
demonstrates  how  the  softening  response  is  affected  by  the  parameter  B,  holding  A  and  Tq  constant. 


As  previously  mentioned,  we  first  fit  our  model  to  the  data  without  including  damage,  i.e.  an  elastic- 
perfectly  plastic  fit  to  TXC  stress-strain  curves.  Then  we  included  damage  in  the  fit  by  iteratively 
adjusting  the  damage  parameters  and  shear  failure  surface  location  to  account  for  pre-peak  plasticity 
and  damage  accumulation.  However,  no  iterative  adjustments  were  made  to  other  aspects  of  the 
fit,  such  as  the  hardening  rule  or  cap  shape  function.  This  is  because  we  wanted  perform  the  Group 
2  lab  test  simulations  first,  to  assess  our  selection  of  the  damage  model  parameters.  The  effect  of 
damage  on  our  original  fits  is  discussed  in  the  following  paragraphs. 

Our  fit  of  the  model,  including  damage,  to  isotropic  compression  and  uniaxial  strain  data  is  shown 
in  Figures  D.IO  and  D.ll.  The  pressure-volumetric  strain  curves  calculated  with  damage  lie  below 
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Strain  (Percent) 

Figure  D.9.  Two  example  fits  with  damage  to  unconfined  compression  data. 

those  of  the  measured  data,  whereas  the  stress  path  with  damage  lies  above  the  measured  data. 
A  better  fit  to  the  data  could  be  obtained  by  iteratively  adjusting  the  hardening  rule,  cap  shape, 
and  damage  model  parameters.  Of  course,  any  adjustment  to  these  parameters  would  also  affect 
our  fit  to  the  TXC  stress-strain  data  previously  shown  in  Figure  D.5.  Thus  an  analyst  must  iterate 
among  numerous  parameters  to  obtain  a  good  fit  of  the  smooth-cap  model  with  damage  to  data. 
One  recommendation  is  development  of  an  automated  optimization  fitting  procedure,  such  as  a 
least-squares  procedure. 

One  final  point  to  note  is  that  our  model  predicts  plastic  volume  expansion,  or  dilation,  which  is 
commonly  observed  in  concrete  test  specimens.  Figure  D.12  demonstrates  that  substantial  volume 
expansion  is  predicted  at  low  confining  pressures.  The  volumetric  strain  in  this  figure  is  negative 
in  expansion  and  positive  in  compression.  Slight  volume  expansion  occurs  at  10  MPa,  but  not  at 
higher  confining  pressures.  The  volume  expands  due  to  large  lateral  tensile  strains  (plastic  flow)  in 
the  TXC  calculations.  Measured  post-pealc  volume  strains  are  not  available  for  comparison  with  the 
calculations  because  they  were  not  provided  by  RDA  due  to  unreliable  strain  gage  measurements. 
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Figure  D.IO.  Fits  of  the  model,  with  damage,  to  pressure- volumetric  strain  data 
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Figure  D.ll.  Fits  of  the  model,  with  and  without  damage,  to  uniaxial  strain  data. 
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Figure  D.12.  Predicted  volume  expansion  under  TXC  loading 
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